Accurate numerical simulations of inspiralling binary neutron stars 
and their comparison with effective-one-body analytical models 
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Binary neutron-star systems represent one of the most promising sources of gravitational waves. In order 
to be able to extract important information, notably about the equation of state of matter at nuclear density, it 
is necessary to have in hands an accurate analytical model of the expected waveforms. Following our recent 
work [1], we here analyze more in detail two general-relativistic simulations spanning about 20 gravitational- 
wave cycles of the inspiral of equal-mass binary neutron stars with different compactnesses, and compare them 
with a tidal extension of the effective-one-body (EOB) analytical model. The latter tidally extended BOB model 
is analytically complete up to the 1.5 post-Newtonian level, and contains an analytically undetermined parameter 
representing a higher-order amplification of tidal effects. We find that, by calibrating this single parameter, the 
EOB model can reproduce, within the numerical error, the two numerical waveforms essentially up to the 
merger. By contrast, analytical models (either EOB or Taylor-T4) that do not incoiporate such a higher-order 
amplification of tidal effects, build a dephasing with respect to the numerical waveforms of several radians. 
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PACS numbers: 04.25.dk, 04.25.Nx, 04.30.Db, 04.40.Dg, 95.30.Sf, 97.60.Jd 



I. INTRODUCTION 

Binary neutron-star inspirals are among the most promising 
and certain target sources for the advanced versions of the cur- 
rently operating ground-based gravitational-wave (GW) de- 
tectors LIGO/Virgo/GEO. These detectors will be maximally 
sensitive during the inspiral part of the signal (around a GW 
frequency of 100 Hz, i.e., significantly below the typical GW 
frequencies at merger, which are around 1000 Hz). The in- 
spiral part of the signal will be influenced by tidal interaction 
between the two neutron stars (NSs), which, in turn, encodes 
important information about the equation of state (EOS) of 
matter at nuclear densities. In other words, the detection of 
GWs emitted from inspiralling NS in the LIGO/Virgo band- 
width could enable us to acquire important information about 
the EOS of NS matter. However, besides getting sufficiently 
accurate GW data from advanced detectors, two conditions 
must be fulfilled for the success of this program: (i) obtain- 
ing a large enough sample of accurate numerical simulations 
of inspiralling binary neutron stars (BNS); (ii) possessing a 
sufficiently accurate analytical model of inspiralling BNS, al- 
lowing the extrapolation of the finite set of numerical sim- 
ulations to the multi-parameter space of possible GW tem- 
plates. Extending the work recently reported in [1], we here 
address issues and provide useful progress on both of them. In 
essence, we will present the results of general-relativistic sim- 
ulations spanning about 20 gravitational-wave cycles of the 
inspiral of equal-mass BNSs, and show how a suitably cali- 
brated effective-one-body (EOB) analytical model of tidally 
interacting BNS systems enables us to accurately reproduce 
the numerically simulated inspiral waveform. 

Numerical simulations of merging BNSs in full general rel- 
ativity have a long history (see the Introduction of [2] for 



a brief review), and the first merger to a hypermassive neu- 
tron star (HMNS) was computed more than ten years ago [3]. 
However, it is only in recent years and with the use of more 
advanced and accurate numerical algorithms that it has been 
possible to obtain a more precise and robust description of this 
process and to include additional physical ingredients such as 
magnetic fields and realistic EOSs. In particular the use of 
adaptive mesh refinement techniques [2, 4, 5] made it possible 
to use very high resolutions, increasing not only the level of 
accuracy, but giving the possibility, for example, to compute 
the full evolution of the HMNS up to black-hole (BH) forma- 
tion [2] or to investigate in detail the development of hydro- 
dynamical instabilities at the time of the merger [2]. Also the 
numerical convergence properties of BNS simulations have 
been studied only recently [6], providing for the first time evi- 
dence of the level of accuracy that it is now possible to achieve 
in the generation of GW templates from these sources. Sev- 
eral groups are now able to simulate BNSs using more realistic 
EOSs (see, e.g., [7-9] and references therein) and to assess the 
possibility to measure their effects in the GW signals. In the 
last two years three different groups were also able to perform 
for the first time the simulations of magnetized BNSs [10-12]. 
One conclusion already reached is that no effect of the mag- 
netic field can be measured in the inspiral waveforms [12], 
while the role of the magnetic field in the post-merger phase 
has been recently investigated in [13] as well as its role in 
the emission of relativistic jets after the collapse to BH [14]. 
Because of their possible connection with the production of 
short gamma-ray bursts (GRBs), numerical simulations have 
also investigated in detail the formation of massive tori and 
their dependence on the initial mass and mass ratio of the bi- 
nary (see e.g., [15]) as well as on the EOS used (see [8, 9] and 
references therein). 
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On the other hand, the program of developing an analyti- 
cal description within general relativity of tidally-interacting 
binary systems has been initiated only recently [16-22]. Over- 
all, this work has brought to light two surprising results. 
First, that the dimensionless expression kg (Love number) 
in the (gravito-electric) tidal polarizability parameter G/if = 
2fc£i?^^+^/(2£ — 1)!! measuring the relativistic coupling (of 
multipolar order £) between a NS of radius R and the external 
gravitational field in which it is embedded strongly decreases 
with the compactness parameter C = GM/{c?R) of the 
NS [18, 19]'. Second, a recent comparison between a numer- 
ical computation of the binding energy of quasi-equilibrium 
circular sequences of BNS systems [23] and the BOB descrip- 
tion of tidal effects [21] suggests that high-order (beyond the 
first order) post-Newtonian (PN) corrections to tidal effects 
tend to significantly increase (typically by a factor of order 
two) the effective tidal polarizability of NSs. 

The main aim of this paper is to present a detailed compar- 
ison between waveforms computed from the tidal-completed 
BOB analytical model of Ref. [21] and waveforms from BNS 
simulations comprising between ~ 20 and 22 GW cycles 
of inspiral [1]. More specifically, we will follow Ref. [21], 
which has proposed a new way of analytically describing 
the dynamics of tidally interacting BNSs, whose validity is 
not a priori limited (like the purely PN-based descriptions 
used in, e.g., [16]) to the low-frequency part of the GW sig- 
nal, but may be extended to higher frequencies, essentially 
up to the merger The proposal of Ref. [21] consists in ex- 
tending the BOB method [24—26], which has recently shown 
its ability to accurately describe the GW waveforms emit- 
ted by inspiralling, merging, and ringing binary black holes 
(BBHs) [27, 28], by incorporating tidal effects in it. We will 
improve the tidally-extended BOB model of Ref. [21] (which 
already contained the IPN contributions to the dynamics) by 
incorporating the IPN contributions to the waveform (from 
[29]), as well as the waveform tail effects (from [30, 31]). 

The paper is organized as follows. In Sec. II we present 
in detail our numerical simulations, briefly reviewing our 
numerical setup, discussing the dynamics of the binaries, 
and presenting the main features of the waveforms. Sec- 
tion III deals instead with the analytical models of the bi- 
nary dynamics and of waveforms that include tidal interac- 
tion (either PN-based or BOB-based). Sec. IV introduces 
some tools, notably a certain intrinsic representation of the 
time evolution of the GW frequency, which is useful for 
doing the numerical-relativity/analytical-relativity (NR/AR) 
comparison. Section V discusses the various errors that af- 
fect the NR phasing. The NR/AR comparison is carried out 
in Sec. VI. We finally present a summary of our findings in 
Sec. VII. Two appendices give additional technical details on 
the use of the waveforms from the numerical-relativity simu- 
lations. 



As a consequence, for a given EOS, the Love numbers of a typical (C ~ 
0.15) NS are found to be about 4 time smaller than their corresponding 
Newtonian estimates, that assume C — >■ 0. 



We use a spacelike signature (— , +, +, +) and (unless ex- 
plicitly said otherwise) a system of units in which c — G = 
AIq ~ 1. Greek indices are taken to run from to 3, Latin 
indices from 1 to 3. 



II. NUMERICAL-RELATIVITY SIMULATIONS 
A. Numerical setup 

The numerical simulations were performed with the set of 
codes Cactus-Carpet-Whisky [32-36]. The reader is re- 
ferred to these references for the description of the details of 
the implementations and of the tests of the codes. Since in this 
work we use the same gauges and numerical methods already 
applied and explained in [2, 6], we also refer the reader to 
these articles for more detailed explanations of the setup only 
briefly recalled below. 

In essence, we evolve a conformal-traceless "3 + 1" formu- 
lation of the Binstein equations in which the spacetime is de- 
composed into three-dimensional spacelike slices, described 
by a metric 7^, its embedding in the full spacetime, speci- 
fied by the extrinsic curvature Kij, and the gauge functions a 
(lapse) and /3' (shift), which specify a coordinate frame (see 
Ref. [34] for details on the latest implementation of the Bin- 
stein equations in the code). For the evolution of the mat- 
ter, the Whisky code implements the flux-conservative for- 
mulation of the general-relativistic hydrodynamics equations 
proposed by the Valencia group [37]. Among its important 
features is that the set of conservation equations for the stress- 
energy tensor T^'' and for the matter current density are 
written in hyperbolic, first-order, and flux-conservative form 
(see Ref. [2] for details on the latest implementation of the 
hydrodynamics equations in the code). 

As initial data we use quasi-equilibrium binaries generated 
with the multi-domain spectral-method code LORENE devel- 
oped at the Observatoire de Paris-Meudon [38]. For more in- 
formation on the code and its methods, the reader is referred to 
the LORENE web pages [39]. In particular, we use irrotational 
configurations, defined as having vanishing vorticity and ob- 
tained under the additional assumption of a conformally flat 
spacetime metric [38]. The BOS assumed for the initial data 
is in all cases the polytropic BOS 

p^Kp^, (1) 

where p and p are the pressure and the rest-mass (baryonic- 
mass) density, respectively. The chosen adiabatic index is 
r = 2, while the polytropic constant is K ~ 123.6 (in units 
where c = G = Mq = 1). For this particular BOS, 
the allowed maximum baryonic mass for an individual sta- 
ble NS is 2.00 A/q, thus leading to a maximum compactness 
A/^Q^,/i? ~ 0.25. The initial coordinate separation of the 
stellar centers in all cases is d = 60 km. 

The physical properties of the two binaries considered here 
are summarized in Table I, where we have adopted the fol- 
lowing naming convention: M%C#, with % being replaced by 
the rounded total baryonic mass AI^^^ of the binary NS sys- 
tem and # by the compactness. As an example, M2 . 9C . 12 is 
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TABLE I. Properties of the binary NS initial data. From left to right the columns show: the name of the model, the total baryonic mass M^ot 
of the system, the total (initial) Amowitt-Deser-Misner (ADM) mass Af^^^j of the system, the total (initial) angular momentum J, the initial 
orbital frequency /orb, the initial maximum rest-mass density /Omax, the mean radius f of each star, the axis ratio A of each star, the individual 
ADM mass M°° of each star as considered in isolation at infinity, the compactness C°° = M'^g /Rns of each star as considered in isolation 
at infinity, the corresponding (quadrupolar) dimensionless Love number k2 and tidal constant K2 as defined in Ref. [21] (see also Eq. (13) 
below). The mean radius is defined as f = {r\~ + + r± + rpoi)/4, where and r-\ are the (coordinate) radii of the star parallel to the line 
connecting the stars, r± is the radius in the equatorial plane perpendicular to that line, and rpoi is the radius perpendicular to the equatorial 
plane. The axis ratio is defined as the ratio between the mean radius parallel to the line connecting the stars and the mean radius in the plane 
perpendicular to that line, namely A = {r± + rpoi)/(rh + r-\). The values of /orb, r. A, M°°, and C°° are computed with the LORENE 
code, the values of Mtot, Af^^^, J, and pmax are instead measured on the Cartesian grid by the Whisky code, and those of k2 (and K2) are 
computed according to Ref. [18]. 



Model 


i\ ,rbar 
*^tot 

(Mq) 


A^ADM 
(Mq) 


J/1049 
(g cm^/s) 


/orb 

(Hz) 


Pmax/lOl-i 

(g/cm^) 


r 

(km) 


A 


M°° 
(M©) 


C°° 


fc2 


T 
Kg 


M2 . 9C . 12 

M3 . 2C . 14 


2.8899 
3.2504 


2.6925 
2.9966 


7.1747 
8.5558 


188.52 
197.03 


4.60 
5.93 


14.2 
13.2 


0.97 
0.97 


1.359 
1.514 


0.1196 
0.1399 


0.09719 
0.07894 


496.09 
183.81 



the binary with total baryonic mass M^^^^ = 2.8899 Mq and 
compactness C = 0.1196. We note that at least as far as the 
tidal effects are concerned, the most important difference in 
the two sets of initial data is represented by the compactness, 
which is smaller in the binary M2 . 9C . 12 than in the binary 
M3 . 2C . 14. Note that the dimensionless EOB parameter 
measuring the strength of the (conservative) quadrupolar in- 
teraction is nearly three times larger for C = 0.12, than for 
C = 0.14. 

The initial data is then evolved either using the (isentropic) 
poly tropic EOS (1) or using the (non-isentropic) "ideal- fluid" 
EOS defined by the condition 

p = pe{T-l), (2) 

where e is the specific internal energy and e = p(l + e) is 
the total energy density. Although these EOSs are idealized, 
they provide a reasonable approximation of the dynamics of 
NSs during the inspiral, so that we expect that the use of re- 
alistic EOSs (with similar compactnesses) would not change 
the main qualitative conclusions of this work. A detailed dis- 
cussion of the consequences of using either EOS will be pre- 
sented in Sec. V. 

As mentioned above, the use of adaptive mesh-refinement 
techniques allows us to reach a considerable level of precision 
and for this we use the Carpet code [33] that implements 
a vertex-centered adaptive-mesh-refinement scheme adopting 
nested grids with a 2 : 1 refinement factor for successive grid 
levels. We center the highest resolution level around the peak 
in the rest-mass density of each star. This represents our rather 
basic form of adaptive-mesh refinement. The timestep on each 
grid is set by the Courant condition (expressed in terms of 
the speed of light) and so by the spatial grid resolution for 
that level; the typical Courant coefficient is set to be 0.35. 
The time evolution is carried out using fourth-order accurate 
Runge-Kutta integration steps. Boundary data for finer grids 
are calculated with spatial prolongation operators employing 
fifth-order polynomials and with prolongation in time employ- 
ing second-order polynomials. 



In the results presented below we have used 6 levels of 
mesh refinement with the finest grid resolution of Amin = 
0.12 A/0 ~ 0.177 km and the coarsest (or wave-zone) grid 
resolution of A„iax = 3.84 Mq — 5.67 km. Each star is com- 
pletely covered by the finest grid, so that the high-density re- 
gions of the stars are tracked with the highest resolution avail- 
able. The refined grids are then moved by tracking the po- 
sition of the maximum of the rest-mass density as the stars 
orbit, and are finally merged when they overlap. In addi- 
tion, a set of refined but fixed grids is set up at the cen- 
ter of the computational domain so as to capture the details 
of the Kelvin-Helmholtz instability (cf. [2]). The finest of 
these grids extends to r = 7.5 A/0 = 11km = 5.52A/ 
for model M2 . 9C . 12 and = 4.95A/ for model M3 . 2C . 14 
(here and in the following A/ denotes the gravitational mass 
of the system at infinite separation, namely the sum of the 
gravitational masses of each NS as computed individually in 
isolation, i.e.. A/ = 2M^g in the notation of Table I). A 
single grid-resolution covers then the region between r = 
I5OA/0 = 221.5km and r = 514.56A/0 = 755.24km 
(or r = 378.63A/ for M2 . 9C.12 and r = 339.87A/ for 
M3 . 2C . 1 4), in which our wave extraction is carried out. The 
resolution is here A = 3.84 A/© = 5.67 km and thus more 
than sufficient to accurately resolve the gravitational wave- 
forms that have initially a wavelength of about 720 km. 

A reflection symmetry condition across the z = plane 
and a 7r-symmetry condition- across the .t = plane are used. 
A number of tests have been performed to ensure that both 
the hierarchy of the refinement levels described above and 
the resolutions used yield results that are numerically consis- 
tent although not always in a convergent regime at the time of 
merger (see the detailed discussion in Ref. [6]). 



^ Stated differently, we evolve only the region {a; > 0, z > 0} applying 
a 180-degree rotational-symmetiy boundary condition across the plane at 
X = 0. 
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B. Overall matter dynamics and gravitational waveforms 

We next briefly recall the physical properties of BNS in- 
spiral and merger as discussed in Refs. [2, 6]. The inspiral 
proceeds at higher and higher frequencies until the time of the 
merger, just before which the stars decompress because of the 
tidal force. At the time of the merger, a Kelvin-Helmholtz in- 
stability develops in the shearing layer formed by the colliding 
stars, which could lead to an exponential growth of magnetic 
fields if these are present [40, 41]; such a large growth was not 
found in recent related works [12, 13], and no magnetic fields 
are included in the simulations reported here. If the total mass 
of the system is sufficiently large, the merged object imme- 
diately collapses to a Kerr BH, while, for smaller masses the 
merger remnant is a HMNS in a metastable equilibrium. Be- 
cause of the excess angular momentum, the HMNS is also 
subject to a bar deformation, being responsible for a copi- 
ous emission of gravitational radiation with peak amplitudes 
that are comparable or even larger than those at the merger 
(cf. Ref. [2]). As the bar-deformed HMNS loses energy and 
angular momentum via GWs, it contracts and spins up, thus 
further increasing the losses. The process terminates when the 
threshold to the collapse to BH is crossed and the HMNS then 
rapidly produces a rotating BH surrounded by a torus of hot 
and high-density material. Although this post-merger evolu- 
tion of the binary is of great interest and is likely to yield a 
wealth of physical information, it will not be further consid- 
ered in the present work, which is instead focussed on the 
analytical modelling of the inspiral phase, up to the merger. 

The GW signal is extracted at different surfaces of constant 
coordinate radius robs by means of two distinct methods. The 
first one is based on the measurements of the non-spherical 
gauge-invariant perturbations of a Schwarzschild BH [42, 43]. 
The second and independent one uses instead the Newman- 
Penrose formalism so that the GW (metric) polarization am- 
plitudes h+ and hx are then related to V'4 by (see Sec. IV of 
Ref. [2] for details of the Newman-Penrose scalar extraction 
in our setup) 

oo e 

h+ - i/ix = V^4 = 51 ^4™ -2Yi,n{0, 0), (3) 

where we have introduced the (multipolar) expansion of ?/'4 in 
spin-weighted spherical harmonics [44] of spin-weight s = 
—2. The coordinate extraction radius is robs = 500 Mq 
for both models, which corresponds to r^hs/M ~ 184.3 for 
M2 . 9C. 12 and to rohs/M = 165.1 forM3 . 2C. 14. The top 
panels of Fig. 1 summarizes most of the information related to 
the I = m = 2 curvature waveforms for the M2 . 9C . 12 
model (left panels) and for the M3 . 2C . 14 model (right pan- 
els). The top panels of the figures show together the modulus 
and the real part of the waveform; the bottom ones, illustrate 
the behavior of the instantaneous GW (curvature) frequency 
AIijj22- Note that the inspiral waveform of M2 . 9C . 12 con- 
tains about 22 GW cycles, while that of M3 . 2 C . 1 4 contains 
about 20 GW cycles. To fix conventions, let us recall that we 
write the waveform as a complex number according to 

V-f" = IVMe-'^'- , (4) 



so that the instantaneous (curvature) GW frequency is sim- 
ply defined as = 0^,„. After the initial junk radiation 
(cf. Ref. [45]) that is responsible for a spike in the modulus 
around t = 200M together with incoherent oscillations in the 
frequency, the complex ■4>'^ waveform becomes circularly po- 
larized (as expected for circularized inspiral), with a modulus 
that grows monotonically in time up to the merger (see upper 
panels of Fig. 1). 

The matter dynamics is reflected in the behavior of the fre- 
quency; for both models we clearly see that LO22 grows mono- 
tonically during the inspiral phase, until it reaches a maxi- 
mum around the "merger". In this work, we phenomeno- 
logically define the "NR merger" as the instant when the 
modulus of the metric waveform /122 (see below) reaches its 
(first) maximum. Roughly speaking, in our simulations the 
"dynamical range" of the dimensionless GW frequency pa- 
rameter il/cL>22 during the inspiral (i.e., before the merger) 
is 0.015 < Muj22 ^ 0.15. Note that, if we were consider- 
ing a conventional 1.4 Mq — 1.4 AIq BNS system, we would 
then have the correspondence /qw/IOOHz ~ 115.4Mw22 so 
that MLd22 = 0.015 corresponds to /gw ~ 173.1 Hz, while 
M1JJ22 = 0.15 corresponds to /gw ~ 1731 Hz. 

In order to perform direct comparisons with (resummed) 
analytical waveforms and since the resummations used in the 
FOB method have been developed (and tested) mainly for 
metric waveforms, we derived the metric waveform by a (dou- 
ble) time-integration of the i/jI^ waveform (the so-obtained 
metric waveform was found to be more accurate than the out- 
put of the gauge-invariant perturbation scheme). We recall 
that the metric waveform is also expanded in spin-weighted 
spherical harmonics with the following convention 

/l+ - i/lx = ^ ^ hiira -2Yhn{d,(l>), (5) 
e=2 m=-l 

SO that the metric multipoles hi„i at time t can be obtained 
from -04™ by double time-integration as 

hi^{t) = f dt' f dt"i:i^{t"). (6) 

J —oo J —oo 

This expression assumes that one knows the curvature wave- 
form on the infinite time interval (—00, t]. Since, however, the 
simulated curvature waveform does not start at an infinite time 
in the past, but at a finite (conventional) time t = 0, one has 
to find a way of determining two (complex) integration con- 
stants accounting from the GW emission from infinite time to 
our present starting time. To do so, we use here an improved 
version of the fit procedure of Ref. [46], which is presented 
in detail in Appendix A. Figure 2 shows the result of this 
process, with the left panels referring to model M2 . 9C . 12 
and the right ones to model M3 . 2C . 14. Note that the wave- 
forms displayed in these figures are obtained from simulations 
with: (i) the non-isentropic (ideal fluid) FOS; (ii) the highest 
available resolution; and (iii) an extraction radius of 500 Mq. 
These will be taken as our fiducial "target" waveforms for our 
NR/AR comparisons, and we will refer to them in the follow- 
ing with the label IFhr500. The numerical uncertainty on 
these target waveforms will be estimated in Sec. V below. 
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FIG. 1. Curvature gravitational waveforms rtp4 (upper panels) and their instantaneous frequency Muj22 (lower panels) for the M2 . 9C . 12 
(left) and M3 . 2C . 14 (right) models. In both cases, the observer's (coordinate) extraction radius is r-Qbs = 500 Mq; this coiTesponds to 
robs/M = 184.3 for M2 . 9C . 12 and Tobs/M = 165.1 for M3 . 2C . 14. 




FIG. 2. Metric gravitational waveforms r/122 and frequencies (upper panels) and the corresponding istantaneous frequency MU22 (lower pan- 
els) obtained from the (double) time-integration of the curvature waveforms of Fig. 1 [see Eq. (6)]. The left panels refer to model M2 . 9C . 12, 
the right panels to model M3 . 2C . 14. The fact that the waveform modulus grows monotonically without evident spurious oscillations is the 
indication of the reliability of the determination of the integration constants. See text for details. 



III. ANALYTICAL MODELS summed EOB description of the conservative dynamics, (ii) 

the resummed EOB description of the waveform, and (iii) one 
of the non-resummed (i.e., PN expanded) descriptions of the 
We recall below some basic information relative to the phasing. 
EOB-based and PN-based descriptions of the binary dynam- 
ics and waveforms that include tidal effects. We follow here 
the general discussion of Ref. [21], to which we refer the 
reader for more details. We consider successively: (i) the re- 
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A. Effective-one-body description of the conservative 
dynamics 

The EOB formalism [24—26] replaces the PN-expanded 
two-body interaction Lagrangian (or Hamiltonian) by a re- 
summed Hamiltonian, of a specific form, which depends only 
on the relative position and momentum of the binary system 
(q,p). For a non spinning BBH system, it has been shown 
that its dynamics, up to the 3PN level, can be described by the 
following EOB Hamiltonian (in polar coordinates, within the 
plane of the motion): 



HEOBir,pr,,p^) = McVl + 2iy{H,ff - 1) , (7) 



where 



H,s=^pl+A{r)i^l + ^ + Z3'-^j. (8) 

Here M = Ma + Mb is the total mass, i' = Ma Mb /{Ma + 
MbY is the symmetric mass ratio, and Z3 = 2i/(4 — 'iv). In 
addition, we are using rescaled dimensionless (effective) vari- 
ables, namely r = rAB(?/GM and = P^c/{GMaMb), 
and pr, is canonically conjugated to a "tortoise" modification 
of r [47]. 

A remarkable feature of the EOB formalism is that the 
complicated, original 3PN Hamiltonian (which contains many 
corrections to the basic Newtonian Hamiltonian i — l/r) 
can be replaced by the simple structure (7)-(8), whose two 
crucial ingredients are: (i) a "double square-root" structure 

Heob ^ 1 + \/p^ + • • • and (ii) the "condensation" of 
most of the nonlinear relativistic gravitational interactions in 
one function of the (EOB) radial variable: the basic "radial 
potential" A{r). The structure of the function A{r) is rather 
simple at 3PN, being given by 



(9) 



where 04 = 94/3 - (41/32)7r2, and u = l/r = 
GM/ (c^rAB)- It was recently found that an excellent de- 
scription of the dynamics of BBH systems is obtained [27] 
by: (i) augmenting the presently computed terms in the PN 
expansion (9) by additional 4PN and 5PN terms; (//) Pade- 
resumming the corresponding 5PN "Taylor" expansion of the 
A function. In other words, the BBH (or "point mass") dy- 
namics is well described by a function of the form 

A°(r) = Pg [1 - 2u + 2vu-^ + a^vu^ + ar-,viC' + a^vu^] , 

(10) 

where P"^ denotes an (n, m) Pade approximant. It was 
found in Ref. [27] that a good agreement between EOB and 
numerical-relativity BBH waveforms is obtained in an ex- 
tended "banana-like" region in the (05,06) plane approxi- 
mately spanning the interval between the points (05,06) = 
(0, -20) and (05,06) = (-36, +520). In this work we will 
select the values 05 = — 6.37, 06 = +50, which lie within this 
region (we have checked that the use of other values within the 
"good BBH fit" region would have no measurable influence on 
our discussion below). 



The proposal of Ref. [21] for including dynamical tidal ef- 
fects in the conservative part of the dynamics consists in sim- 
ply using Eqs. (7)-(8) with the following tidally-augmented 
radial potential 



A{u) = A'\u) + A^^'^'^^u) 



(11) 



Here A^{u) is the point-mass potential defined in Eq. (10), 
while ^'idai j-y-j ^ supplementary "tidal contribution" of the 
form 



A 



tidal 



l>2 



(12) 



where the terms kJv?'^'^'^ represent the leading-order (LO) 
tidal interaction, i.e., the Newtonian order tidal interaction. 
The dynamical EOB tidal coefficients kJ are functions of the 
two masses Ma and Mb, of the two compactnesses Ca.b = 
GMa,b/ Ra,b, and of the two (relativistic) Love numbers 



r A,B 



of the two objects [18-20]: 



T 
Kg 



Mb Mf 



\Ma + MbY'+^ Cf+ 



j + {a ^ b} 



1 



22£-l Q2e+1 ' 



(13) 



where the second line refers to an equal-mass binary, as the 
ones considered here. Note in Table I the rather large numer- 
ical values for the £ = 2 tidal coefficients: k^(C = 0.12) ~ 
496 and k'^{C = 0.14) ~ 184. In our EOB modefling we also 
use the higher multipolar tidal coefficients and kJ, which 
are even larger than (e.g., kJ{C = 0.12) ~ 20318), al- 
though their effect is marginal in view of the higher power of 
u (namely u^''+^) with which they enter the A{r) potential. 

The additional factor A^"^'^'(u) in Eq. (12) represents the 
effect of higher-order relativistic contributions to the dynam- 
ical tidal interactions: next-to-leading-order (NLO) contri- 
butions, next-to-next-to-leading-order (NNLO) contributions, 
etc. Here we will consider a "Taylor-expanded" expression 



A 



tidal 



(u) 



-(f) 

a\ u 



+ a, u 



(14) 



where a-lf'' are functions of Ma, Ca, and k^ for a general bi- 

nary. The analytical value of the (£ = 2) IPN coefficient a\ 
has been reported in [21] (and recently confirmed in [48]). In 
the equal-mass case, it yields a^'' = 1.25. By contrast, there 

are no analytical calculations available for a\ with £ > 2, 

-it) 

nor for the 2PN tidal coefficients 0.2 ■ Indeed, one of the 
main aims of the present work will be to constrain the value 
of cij^'' by comparing the EOB predictions to numerical data. 



(2) 



Effective-one-body description of the waveform and 
radiation reaction 



Let us first recall that the EOB formalism defines the radia- 
tion reaction from the angular-momentum flux computed from 
the waveform. Concerning the waveform, in the case of BBH 



7 



systems, the EOB fonnalism replaces the PN-expanded mul- 
tipolar (metric) waveform /i™ by a specifically resummed 
"factorized waveform" [31, 49], say h^^^ (where the super- 
script is added to signal the absence of tidal effects). This 
tidal-free multipolar waveform includes resummed ver- 
sions of very high-order PN effects in the phase and the modu- 
lus, in particular tail effects. Actually, in the present work, we 
have used a factorized waveform which includes in the modu- 
lus (but not in the phase) the new (5PN accurate) ly ~ terms 
recently computed in [50]"'. We also included in the two 
next-to-quasi-circular parameters (01,02) as in Ref. [27]"*. 

When considering tidally interacting binary systems, one 
needs to augment the BBH waveform by tidal contribu- 
tions. Similarly to the additive tidal modification (11) of the 
A potential, we will here consider an additive modification of 
the waveform, having the structure 



hlrn — hg 



'tidal 



(15) 



This is slightly different from the factorized form introduced 
in Eq. (71) of [21] and used in [1]. The above additive form 
turns out to be more convenient for incorporating higher-order 
relativistic corrections to the tidal waveform. Using the recent 
computation [29] of the IPN-accurate Blanchet-Damourmass 
quadrupole moment [51] of a tidally interacting binary sys- 
tem (together with the Newtonian-accurate spin quadrupole 
and mass octupole) and transforming their symmetric- trace - 
free tensorial results into our £m-multipolar form, we have 
computed the corresponding IPN-accurate value^ of /i22'^', Taylor 
as well as the OPN-accurate values of /il^''', and °-3-5 i^) ^ ^ [ 



here as an effective parametrization of all the higher-order ef- 
fects not covered by the current analytical knowledge (both in 
the conservative dynamics and in the radiation reaction). Note 
also that, while in the general case such a parameter should be 
allowed to depend on the mass ratio and the compactnesses, in 
the equal-mass case that we consider here, it is a pure number. 
We will use below the comparison between NR simulations 
and EOB predictions to constrain the value of the effective 
higher-order parameter 0:2 . 



C. PN-expanded Taylor-T4 

Tidal effects can be accounted for also via modifications of 
one of the non-resummed PN description of the dynamics of 
inspiralling binaries [7, 16, 20]. Reference [20], in particular, 
has recently suggested to use as baseline a time-domain T4- 
type incorporation of tidal effects. We recall that the phasing 
of the T4 approximant is defined by the following equations 



^^^22^ _ r, 3/2 

dx 64 ( • 
— = —vx \a 



3.5 



\x)+a'^^^\x)} , 



(16) 



where aj^^"'^ is the PN expanded expression describing point- 
mass contributions, given by 



'21 ' "33 

^tidai addition, using the general analysis of tail effects 
in Refs. [30, 52] and the resummation of tails introduced in 
Refs. [31, 53], we were able to further improve the accuracy 
of these waveforms by incorporating (in a resummed manner) 
the effect of tails (to all orders in M). From a PN point of 
view, this means, in particular, that the tidal contribution we 
use to the total metric waveform is 1.5PN accurate. 

In summary, the EOB tidal model that we use here is ana- 
lytically complete at the 1.5 PN level. In addition, we adopt 
the simplifying assumption that the higher-multipolar tidal- 
amplification factors ^^"^'^'(it), for ^ > 2, are taken to coin- 
cide with the £ = 2 one. This means that the EOB model that 
we will use here contains only one (yet undetermined) higher- 
order flexibility parameter, say a2, that is taken to replace the 



743 
336 



— 1^ ] X + A-KX^^'^ 

4 



/ 34103 13661 
V 18144 ^ 2016 ' 
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16447322263 1712 



139708800 
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56198689 
217728 ' 
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tt\ , 856 , , 

-(256 + 451Z.)-— ln(16x) 



541 
896* 

x' 



4032 



358675 
6048 



91495 
1512 



] nx'/^ 



(17) 



and where a*"*''' is the tidal contribution. From [29] the latter 
is given at IPN accuracy by 



various a 



^tidal 



with £ = {2, 3, 4, . . . }, entering Eq. (14), i.e. 



(x-)= J2 ciLo{Xi)x^{l + ai{Xi)x) , (18) 



a2 (and, similarly, a\ 



T,(2) - 



ai). Note that. 



although this parameter is formally of 2PN order, it is used with 



I=A.B 



Ah 



jl2 - llXi 

~x'j 



^ As in Ref. [49] we resum the £ = 2,m = 2 modulus by using the Pade- 
resummed function /|'2^(x; u) = Pj u)]. 
Since both M2 . 9C . 12 and M3 . 2C . 14 are equal-mass binaries, we fix 
ai = —0.0439 and 02 = 1.3077, according to the EOB/NR comparison 
(for a BBH equal-mass system) of Ref. [27] . 

^ We leave a detailed presentation of our results to future work. Let us 
however mention that, notwithstanding some statements in footnote 4 of 
[29], the IPN-accurate (circular) quadrupolar waveform exactly matches 
the form given in Eq. (71) of [21] (which was expressed in terms of 
frequency-related gauge-invariant quantities). 



and 



ai{X) = 



4421 - 12263X + 26502^2 - 18508^3 



336(12- IIX) 
where we introduced the auxiliary quantity 



— \ 



I ^A,B. 



(19) 



(20) 



(21) 
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In the particular case of equal-mass binaries, Xa = Xb 



X = 1/2, Ca 
has the form 



Cf 



^tidal / 



C, and the tidal contribution 



{x) =26k^x'^ (l + a^'^x) 



-^tidal 



(22) 



with aj^ = 5203/4368 « 1.19. 

Similarly to the inclusion of yet undetermined higher-order 
effects in the tidally-augmented EOB formalism via the effec- 
tive parameter 0.2, we will consider below an effective modi- 
fication of the IPN result (22) of the form 



^tidal 



(x) 



26 K^x^ (1 



T4 

+ a, X + a. 



T4^2 



x^), (23) 



with an effective higher-order parameter^ aj^, which we will 
constrain by comparing NR data to the T4-predicted phasing. 

Let us mention that, in the case of inspiralling BBH sys- 
tems, several studies [31, 46, 54] have shown that the nonre- 
summed Taylor- T4 description of the GW phasing was signif- 
icantly less accurate than the EOB description, especially for 
mass ratios different from one. Ref. [21] has also shown that, 
in the presence of tidal effects, it was predicting GW phases 
that differed by more than a radian with respect to the tidal- 
completed EOB model. Below, we will investigate how the 
T4 phasing based on Eq. (16) differs from the EOB one, both 
in the absence (Eq. (22)) and in the presence (Eq. (23)) of the 
higher-order parameter aj^. 



IV. CHARACTERIZING THE PHASING: THE (oj) 
FUNCTION 

In order to measure the influence of tidal effects it is useful 
to consider the "phase acceleration"^ lj = doj/dt = <:P(p/dt^ 
as a function of u, say lu = a{uj) (here ui = UJ22 can be either 
the curvature or the metric instantaneous GW frequency). In- 
deed, as emphasized in [31], the function a{uj) is independent 
of the two "shift ambiguities" that affect the GW phase (f>{t), 
namely the shifts in time and phase. The a{Ld) diagnostics (es- 
pecially in its Newton-reduced form = a{uj)/{c^uj^^/'^), 
with = ^2^/'^i^, is a useful intrinsic measure of the qual- 
ity of the waveform and it has been used extensively in recent 
analyses of BBHs [46, 53, 55, 56]. 

Here we will use another dimensionless measure of the 
phase acceleration: the function Q^{uj). It is defined as the 
derivative of the (time-domain) phase with respect to the log- 
arithm of the (time-domain) frequency: 



d In 



Lud(j)/dt ur' up' 
duj/dt u! a{u!) 



(24) 



^ We found that the 1.5PN fractional contribution a^^^x^^^ to a*"*^'(x), 
predicted by our 1.5PN-accurate EOB waveform, has (like the IPN contri- 
bution) only a small effect on the phasing compared to the large amplifica- 
tion that we will need to agree with NR data. This is why we only consider 
here, for simplicity and for easier comparison with the 2PN EOB parameter 
02 , the formally 2PN parameter a J'' . 

^ In the text of this Section t and lu denote the dimensionless quantities i = 
t/M andtj = Mu. 
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FIG. 3. Exploring the properties of curves computed within the 
EOB model for three binary systems. Tidal interactions are approx- 
imated at LO. The inset shows a magnification, in order to highlight 
the differences among the curves. 



Note that, as a consequence of this definition, the (time- 
domain) GW phase </>((^i,(^2) accumulated between frequen- 
cies (wi, W2) is given by the following integral: 



Q^dlnu) . 



(25) 



Stated differently, the function (w) measures the number 
of GW cycles spent by the binary system within an octave of 
the GW frequency cj (it is therefore analogous to the "quality 
factor" Q of a damped oscillator). Let us also note that, in the 
stationary phase approximation, Q^^ enters as an amplification 
factor of the signal, so that the squared signal-to-noise ratio is 
equal to [57] 



P 



= 4 / din 



(26) 



where A denotes the ampHtude of the time-domain met- 
ric waveform and where Sn{f) denotes the one-sided noise 
power spectral density and / = uj/{2tt). 

In view of its definition, is a useful quantitative indi- 
cator of the physics driving the variation of uj. Indeed, a 
change of Qw(w) of the order ±1 during a frequency "octave" 
ln(w2 /ui) ~ 1 corresponds to a local dephasing (around w) of 
Afp ~ ±1 rad. Because such a dephasing (if it occurs within 
the sensitivity band of the detector) can be expected to sig- 
nificantly affect the measurability of the signal, it is probably 
necessary to model Q^^ with an absolute accuracy of about ±1 
(see Ref. [55] for a quantitative discussion of the admissible 
error level on in the BBH context). 

We start our analysis by comparing the functions (as 
predicted by the EOB formalism) for the (metric) gravitational 
waveforms /i22 generated by three (equal-mass) binary mod- 
els, namely a BBH and the two BNS systems discussed in 
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FIG. 4. Obtaining the Qu) diagnostic from a suitable fitting procedure of the GW phase (for both curvature and metric waveforms). The two 
vertical lines on the left panels indicate the time interval At/M = [1000, 2290] where we fit the NR phase with Eq. (28). For complete- 
ness we also display the real part of the metric waveform. On the right panels, the (light) dashed lines refer to the Quj obtained by direct 
numerical differentiation of the raw data; the solid lines are instead obtained from the fitted phase. Although the curves displayed here refer to 
model M3 . 2C . 1 4, similar results are obtained also for the binary M2 . 9C . 12. 



Sec. II A. To simplify the discussion, these functions are com- 
puted with the LO tidal interaction A({u) = 1. [We will sep- 
arately study below the effect of changing Ag{u).] 

Figure 3 compares the properties of the functions by 
showing together the curves for the three binaries versus their 
corresponding GW frequency. A number of remarks are worth 
making. First, Q^j is a large number that diverges in the 
small-frequency limit. This follows from the fact that in 
the limit uj ^ one has a{uj) ^ c^w^^/'^, and then, via 
Eq. (24), = l/(c,.w5/3) {c/vf. Second, the pres- 
ence of tidal interactions decreases the "point-mass" value of 
Qu, by an amount that is (essentially) proportional to kJ. In 
other words, tidal effects "accelerate" the inspiral by reducing 
the number of cycles spent around a given frequency. In par- 
ticular, BBHs (which have vanishing tidal constants [18, 19]) 
are effectively the binaries that spend the largest time at any 
given frequency. Finally, note that since Quj is a large num- 
ber, the fact that the curves look relatively close on the large- 
scale plot can be misleading, since the corresponding accumu- 
lated relative phase difference can actually be large (see inset, 
which shows that the absolute differences between the various 
Quj{^^) is of order 10, corresponding to integrated dephasings 
of order 10 radians). 

Although the calculation of the phase "quality-factor" Q^^ is 
straightforward within the BOB framework, this is not the case 
when is to be calculated from the NR (either curvature 
or metric) waveforms. Indeed, the direct computation of the 
functions from raw data is in general made difficult by 



the presence of both high-frequency noise in uj{t) and of low- 
frequency oscillations probably due to a residual eccentricity. 
This is illustrated in the right panels of Fig. 4, where we show 
with (light) dashed lines the raw NR functions obtained 
by direct time-differentiation of the NR curvature (top panel) 
or metric (bottom panel) phase for the binary M3 . 2C . 14. A 
fourth order accurate finite differencing algorithm has been 
used to compute the derivatives. Similar results have been 
obtained also for the binary M2 . 9C . 12. 

The right panel of Fig. 4 shows that the two time-derivatives 
involved in the definition of Quj{(^) amplify considerably the 
high-frequency noise contained in the NR phase evolution, 
and make it impossible to extract a reliable value of Qui^^) 
from such a direct numerical attack. To tackle this prob- 
lem, one needs to filter out the high-frequency numerical 
errors in the time-domain phase before effecting any time- 
differentiation. To do this, we found useful to "clean" the 
phase (j){t) by fitting the NR phase to an analytic expression 
that is inspired by the PN expansion. More precisely, after 
introducing a formal "coalescence" time tc and defining the 
quantity 



-1/8 



(27) 



we fitted the time-domain NR phase <j) [t) to an expression 
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of the form 

2 _r 

0(^;ic,P2,P3,P4,0o) = 00 + --a; 

X {I + P2X^ + P3X^ + pix"^) . (28) 

In this expression, we have set the lower coefficients po and 
pi to Po = 1 and pi = 0, as suggested by the correspond- 
ing lowest-order PN expression (see, e.g., Eq. (234) of [58]), 
but we left tc, (po, and the higher-PN p/s as free coefficients 
to be determined from the NR data. The basic idea is that of 
using a simple analytical form that incorporates the leading 
trend of to remove the influence of the numerical errors 
while leaving some flexibility in the subleading part of the 
phase evolution that is influenced by tidal effects. We view 
the fitting parameters P2,P3,P4 as effective parameters for de- 
scribing tidal-phasing effects. 

Such a fit of the phase evolution can be reliably done only 
in a limited time interval. Indeed, one has to cut off both the 
early phase of the inspiral (where the numerical data is too 
noisy), and the last few cycles before the merger (where the 
PN-based fit is no longer a good approximation). We present 
in Appendix B a detailed discussion of the optimal choice of 
the time interval where to make the fit, as well as a series 
of consistency checks. See also the discussion at the end of 
Sec. VB. 

Let us start by discussing the application of this procedure 
to the GW phase (both curvature and metric) of the binary 
model M3 . 2C . 14. The result of this fitting is shown by the 
solid lines in the right-panels of Fig. 4 (top, curvature phase; 
bottom, metric phase). The time interval on which we could 
rehably apply the fitting procedure is It/M = [1000,2290]. 
This time window is indicated by the dashed vertical lines in 
the top-left panel of Fig. 4, were we show together the time 
evolution of both the curvature (dashed, red online) and met- 
ric (solid) GW frequencies. For completeness, the lower-left 
panel of the same figure translates this information in terms 
of GW cycles of the metric waveform. Note that this time in- 
terval excludes the first 4 GW cycles (whose NR frequency is 
indeed seen to be quite noisy), but covers about 10 GW cy- 
cles, and ends around 2 GW cycles before the merger (defined 
as the maximum of the modulus of the metric waveform; the 
modulus of the metric waveform is indicated by a dashed line 
on the left-bottom panel of the figure). The corresponding 
frequency interval can be visualized on the right panels, and 
is listed in the fifth column of Table 111. Similar results are 
obtained also for the M2 . 9C . 12 data (see Fig. 10 below). In 
this case, the time interval we use is It/M = [1300,3366], 
with the corresponding frequencies listed in the seventh col- 
umn of Table 111. Finally, notice that for this model the inspiral 
is longer than in the previous case and so this interval actually 
corresponds to 14 GWs cycles. In addition, similarly to the 
other case, our choice of fitting interval excludes the first 5.5 
GW cycles, and ends about 2 GW cycles before merger 

As we will discuss below, although the frequency windows 
where our cleaning procedure allowed us to compute an esti- 
mate of the NR (a;) functions do not cover the full inspiral, 
these estimates will give us access to important information 
for performing quantitative comparisons with the predictions 
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FIG. 5. Comparing waveforms from isentropic (dashed) and non- 
isentropic (solid) evolution for BNS model M3.2C.14. Wave- 
forms are computed with the highest resolution and extracted at 
''obs = 500 Mq. The corresponding phase difference (^P°'yHR500 _ 
^ifhrSoo displayed in Fig. 6. 



of the FOB (and Taylor T4) analytical models. 



V. NUMERICAL ERROR-BUDGET 

The aim of this section is to discuss the various errors af- 
fecting the numerical waveforms extracted (for both models) 
at 500 il/0 and computed with the highest resolution. Such 
a discussion will in turn allow us to estimate an uncertainty 
range on the analytical parameter 0.2 representing the not-yet- 
calculated, high-PN-order tidal effects entering the BOB de- 
scription of the phasing. 

We will discuss in turn the numerical errors entailed by 
three different effects: (i) the choice of EOS (isentropic ver- 
sus non-isentropic evolution); (ii) the finite extraction radius; 
(iii) the finite resolution. We will perform this analysis both 
by comparing waveforms in the time domain and by means of 
the diagnostic. 



A. Time-domain analysis 

1. Non-isentropic evolutions 

As discussed in Sec. 11 A, we have evolved the binaries us- 
ing either a (isentropic) polytropic EOS or a (non-isentropic) 
ideal-fluid EOS. We recall that, in the absence of large-scale 
shocks (like those taking place at the merger), the two EOSs 
are equivalent and should therefore yield evolutions that dif- 
fer only at machine precision. In practice, however, when 
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using the ideal-fluid EOS small shocks are produced in the 
very low-density layers of the stars even when these orbit [2]. 
These small shocks channel some of the orbital kinetic en- 
ergy into internal energy, leading to small ejections of matter 
(i.e., amounting to a total of ~ 10^^71/©), and are thus re- 
sponsible for slight differences even during the inspiral. Since 
we are here presenting the results of simulations that are con- 
siderably longer than any presented so far and in particular 
of those in Refs. [2, 6], it is important to quantify the influ- 
ence of these non-isentropic effects. Concentrating on model 
M3 . 2C . 1 4, we show in the top-panel of Fig. 5 the real parts 
of the rV-'p waveforms computed with the two EOSs as ex- 
tracted at robs = SOOA/q = 165. lA/. The bottom panel dis- 
plays the corresponding instantaneous frequencies for com- 
pleteness. As customary in comparing waveforms in the time 
domain, one allows for arbitrary relative time and phase shifts 
(r, a). These quantities can be determined in various ways, 
for example by means of the two-frequency pinching tech- 
nique of Ref. [59]. In this paper we find it useful to use the 
method used in Ref. [54] to compute (r, a). More precisely, 
given two timeseries of the phase {(j)i{ti), (^2(^1)} defined on 
a given time interval [t^ , tn] that is covered by N numerical 
points ti, with i = 1, 2, . . . , A^, we define the quantity 



A(j>{U,T,a) = 02 (ij 



a 



(29) 



and determine t and a such that they minimize the "reduced" 
quantity 



1 ^ 

-^(A</>(i„r,a))= 



(30) 



The minimization on a is done analytically, while that on r is 
done numerically. Note in addition that the square root of the 
minimum value of Eq. (30), say 



\ 



1 ^ 

-^(A0(i„T,a)) 
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min 



(31) 



i=l 



has the meaning of a root-mean-square deviation of the phase 
difference A0 over the interval [tL,t]i]; as such, it can also 
be employed to give a quantitative measure of a phase dif- 
ference (and thereby of some phase errors).*^ The phase dif- 
ference A(j){t) = (t>2{t) - Mt) = (AP°lyHB500 _ 0IFhr5OO 

(least-square minimized on the time interval [tL,tii]/M = 
[300, 2540]) is represented as a dash-dotted line (solid light 
blue) in Fig. 6. One sees that the instantaneous phase differ- 
ence varies roughly between +0.2 rad and —0.1 rad on this 
time interval, which corresponds to a "two-sided" [59] phase 



uncertainty of the order =^-(0.2 - (-0.1)) = ±0.15 

rad. The information of Fig. 6 is completed by Table 11, where 
we list the norm of the phase difference [i.e., the maximum 
absolute value of A0(<)], labelled ||A0||°°, the root-mean- 
square (Ta0 as computed above, and the corresponding time 
interval [tL,tR] that is used to compute (a, r). Note that aA<p 
gives a measure of the phase difference which is always sig- 
nificantly smaller than the £°° norm. Indeed, these two quanti- 
ties measure different aspects of a phase difference, and, when 
the time variation of A0(t) is dominated by low-frequency ef- 
fects (which can be roughly modelled as power laws), the av- 
eraging involved in the definition of 0-^4, will lead to a small- 
ish ratio crA0/||A(/)||°° < 1 linked to integrals of the type 
j^dte'^ = 1/(271+1). 



2. Finite-radius extraction 

We next discuss the phasing error introduced by the fact 
that our high-resolution target waveforms, for both models, 
are extracted at the finite coordinate radius robs = 5001^0. 
Note that, when expressed in units of the gravitational mass 
M of the binary at infinite separation, this value corresponds 
to robs = 134.9A/ for M2 . 9C . 12 and robs = 165.1 A/ for 
M3 . 2C . 14, i.e., , for one model waves are actually extracted 
slightly farther than for the other For both models we have at 
our disposal several extraction radii, so that we can estimate 
the phasing error linked to the finite extraction radius as fol- 
lows: (i) We used the raw data extracted at radii r = 
{400, 450, 500} Mq; (ii) We time-shifted them so that this 
triplet of timeseries is expressed as a function of the (coordi- 
nate) retarded time = i— ?■- 2A/adm In [r/(2A/ADM) — 1]; 
(iii) We separated each curvature waveform in phase and am- 
plitude as functions of u = l/r (cf. page 6); (iv) We fitted 
each resulting triplet of timeseries to a linear polynomial in 
the triplet of inverse extraction radii: c°° (<* ) + ci (<* ) /r. The 
quantities c°°(t*) [i.e., A°°{t^,) and (f>°°{t^,)] yield estimates 
of the amplitude and phase of the infinite-radius extrapola- 
tion of ripj^. We then compare the radius-extrapolated phase 
4>°° (t* ) to the phase extracted at the outermost radius, allow- 
ing for additional time and phase shifts (which are determined 
by the least-square minimization discussed above). 

The time evolution of the phase differences computed in 
this way are shown in Fig. 6 for model M3 . 2 C . 1 4 (top panel, 
dash-line) and for M2 . 9C . 12 (bottom panel). This local in- 
formation is completed by the "global" quantitative informa- 
tion (| I A0I |°°, (TA^) listed in the last two columns of Table 11. 
On the basis of this analysis, we estimate that, for both mod- 
els, the phase uncertainty due to finite extraction is of order 
A(j) w ±0.05 rad almost up to the merger, i.e., roughly 100 M 
before the peak of the GW frequency. 



' We note in passing that the alignment procedure also highhghts the weak 
dependence on the EOS of the late part of the waveform: although the 
inspiral of the non-isentropic waveform is about 150 A/ longer than the 
con'esponding isentropic one, the growth of Mu)22 (and the con'esponding 
phasing) is qualitatively and quantitatively very close for both models until 
Mlo22 peaks for the first time. 



3. Finite-resolution error 

Finite-resolution errors have already been discussed in de- 
tail in our previous work [6], which used the same numeri- 



12 




IF: Radius extrapolated 

IF: Resolution extrapolated @2OOjl/0 

- - Poly @500Mq 
Poly: radius extrapolated 



1000 



t/Af 



1500 



2000 



TABLE II. Uncertainty estimates on tiie piiase (in radians) of ripi^, 
computed in the time domain, for both models. From left to right, the 
columns report: the model name, the EOS, the coordinate extraction 
radius, the type of extrapolation that is performed on the waveform 
(either in extraction radius or resolution), the time interval on which 
the of the phase difference is minimized, the norni of the 
phase difference over this interval, and the root-mean-square of the 
phase difference. 
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FIG. 6. Estimate of the phase uncertainty in the time domain for 
model M3 . 2 C . 1 4 (top) and M2 . 9C . 1 2 (bottom). The figure shows 
the phase difference between different "post-processed" numerical 
curvature waveforms rtjj'p (in particular, extrapolated in resolution 
and/or extraction radius) and the one obtained with the ideal-fluid 
EOS and extracted at robs ~ 500 A/© . 



cal setup (i.e., the same resolution and grid structure) adopted 
here. Skipping the details, we recall that it was shown there 
that, at the resolution that we are using in this work, the dy- 
namics and waveforms are in a convergent regime, with a con- 
vergence rate a that is ~ 1.8 during the inspiral phase and 
drops to ~ 1.2 after the merger, when large-scale shocks ap- 
pear As the computational cost of the calculations presented 
here is already at the limit of what can be reasonably af- 
forded, we have decided to estimate the truncation-error of our 
present waveform by assuming that the inspiral convergence 
rate cr ~ 1.8 found in our previous work [6] approximately 
holds in the present (numerically similar) case. We have then 



selected the more compact binary M3 . 2C . 14 and used only 
two simulations with different resolutions. More specifically, 
we have considered a "high-resolution" simulation, where 
the finest refinement level has a resolution Ah = 0.12 Mq, 
and a "low-resolution" simulation, with Al = 0.15 Mq. 
For this particular comparison the waveforms are extracted 
at Tobs = 200 A/q. When comparing the low- and high- 
resolution curvature waveforms, after suitable (t, a) align- 
ment, one discovers that the phase difference accumulated be- 
tween the two resolutions over a timescale of 2300M during 
the inspiral is about 0.45 rad (corresponding to a relative error 
of ~ 0.36%). Using the convergence rate discussed above, 
we can now Richardson-extrapolate the results obtained with 
the two resolutions and obtain an estimate of the "infinite- 
resolution" waveform. More precisely, we model the suitably 
aligned, low- and high-resolution phase evolutions as 



bAAt) = Mt) + km^nY 



(32) 
(33) 



where (j)o{t) represents the infinite-resolution phase (A 0). 
From the above equations, we obtain the following estimate 
of the infinite-resolution extrapolation of the phase evolution 



Mt) 



(AL)"0AH(<)-(AHr0Ajt) 

(Al)- - (Ah)- 



(34) 



We performed the same extrapolation also on the waveform 
modulus, so as to have access to the complete extrapolated 
curvature waveform. The solid line in Fig. 6 displays the 
phase difference i/iq^^™ — ^iFhrSOO jj^jg indicates a phase 
uncertainty of Acj) w ±0.5 rad on ^i^hrSOO measured up 
to about lOOAf before the maximum of Muj22- See Table 11 
for the corresponding global measures of the phase uncer- 
tainty, ||A(/)||°° and cta^- Note that these uncertainty esti- 
mates are much larger than those normally computed for bi- 
nary BH simulations for the same computational costs (see, 
for instance, [60]). This is the natural consequence of the 
smaller resolution employable here and of the lower-order 
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FIG. 7. Left panel: span of Qtj's due to the various approximations to ttie curvature wavefoims from model M3 . 2C . 14. Right panel: the 
corresponding differences AQ^^ = — Q^J^^'''^^' between the various curves and the fiducial one obtained from the phase computed at the 
highest resolution and extracted at 500 Mq . 



TABLE III. Uncertainty estimates on the r'i/'l phase of the IFhr500 fiducial simulations obtained from integration of the differences between 
Quj's. From left to right the columns report: the model name, the EOS, the coordinate extraction radius, the type of extrapolation that 
is performed on the waveform, the frequency interval MI^ where the cleaning procedure is applied, the corresponding time interval It, the 
accumulated phase difference A(f>^^ = — c/j'^^hrSoo ^ common frequency interval MIZ, the number of GW cycles on the same frequency 

interval, and the relative phase difference A(^^^ = A(^^^/<^^,^. We choose the common interval of integration to be MIZ ~ [0.045, 0.067] 
for model MS . 2C . 14 and MI^ = [0.037, 0.054] for model M2 . 9C . 12. 
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convergence that one achieves when solving the hydrodynam- 
ics equations. Since this error is deduced only after assum- 
ing a certain convergence order (in addition obtained from 
a sUghtly different numerical setup) it will be used below 
only to estimate a rough uncertainty range on the value of 
the higher-order EOB tidal correction parameter a2- We will 
comment more on this in the next Sections. 

Adding in quadrature the various uncertainties computed so 
far to obtain a total error bar on the phases of the IFhr 500 
data for the M3.2C.14 model would give a (two-sided) 
time-domain phase uncertainty A(/) ~ ±\/0.15^ + 0.05'^ ~ 
±0.16 rad, when excluding the uncertainty due to the finite 
resolution, or Acj) ~ iVO.lS^ + O.OS^ + 0.5^ ~ ±0.52 rad 
when including it. Alternatively, if we add in quadrature the 
root-mean-squares of the corresponding phase errors we find 
aA4, — ±0.07 rad, when excluding the uncertainty due to 
the finite resolution, and cta^ — ±0.32 rad when including it. 



Clearly the resolution-extrapolation error is dominating the er- 
ror budget. In view of the uncertainty in estimating this source 
of error, we will not directly use these time-domain phase- 
error levels in estimating the uncertainties in the comparison 
between the EOB, T4, and NR phasings. As we will discuss 
next, we prefer to express the information gathered above on 
numerical errors in terms of the corresponding Q^^ curves. 

B. Qaj analysis 

In Sec. IV we have introduced = w^/w as a convenient, 
intrinsic diagnostics to describe the phasing of the waveform. 
In particular, it allows us to better visualize the influence of 
tidal effects on the phasing, as well as to compute the dephas- 
ing accumulated on a given frequency interval. It is then use- 
ful to recast the various time-domain phase uncertainties on 
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FIG. 8. Sensitivity of Qui to the phase model used in the fitting pro- 
cedure. Note that the n = 4 and n = Q curves are barely distinguish- 
able on the plot. See text for details. 



the high-resolution wavefonn extracted at 500 M© discussed 
above, in terms of Q^j- In practice, we apply the cleaning 
procedure on each waveform of Table 11 so as to obtain four 
Qi^ curves. These curves are displayed together in the left 
panel of Fig. 7, while the fifth column of Table 111 lists the 
specific frequency intervals that were selected to apply the 
cleaning procedure. For a quantitative assessment of the dif- 
ferences between the curves, we present in the right panel 



of Fig. 7 the quantity AQ, 



where the labelling X indicates any curve other than our fidu- 
cial one, IFhr500. Although the information conveyed by 
this figure is qualitatively analogous to the time-domain anal- 
ysis. Fig. 6, it is made here both independent of any phase- 
alignment procedure and simpler to quantify. First of all, the 
figure shows that the extrapolations in radius and in resolu- 
tion act in different directions: the first one pushes the curve 
down (i.e., less GW cycles accumulated on a given frequency 
interval, tidal effects look stronger), while the second one 
pushes the curve up (i.e., more GW cycles accumulated and 
tidal effects look weaker). This result is qualitatively com- 
patible with the corresponding curves in Fig. 6, whose 
slopes have opposite signs. In addition, by integrating in fre- 
quency the AQuj curves on the common frequency interval 
MI^ ~ [0.045,0.067] one obtains an estimate of an actual 
accumulated phase error that can be compared to our previous 
time-domain results (i.e.. Fig. 6). The result of this integra- 
tion is given in the seventh column of Table 111. Note that the 
A0^4 computed in this way is typically significantly larger 
than what was estimated above in the time domain. For in- 



stance, regarding the comparison with the resolution extrap- 
olated waveform, the -based procedure indicates a phase 
difference of about 1.3 rad over by contrast, inspecting 
Fig. 6, where the vertical (red) dashed line corresponds to 
in the time-domain, we read from the plot an accumulated 
phase difference on this interval of about 0.8 rad, i.e., about 
40% smaller Similar results hold for the other phase compar- 
isons. This increase in the estimated phase errors is probably 
due to the additional uncertainty brought by the necessity to 
use a phase-cleaning procedure to compute each Q-^{lo) (see 
below). This is the price we have to pay to be able to have the 
convenience of an intrinsic diagnostic of the phase evolution. 

A separate discussion is needed when comparing isentropic 
and non-isentropic Q^j curves. Figure 7 indicates that the 
curve corresponding to the ideal-fluid EOS lies above the 
polytropic one, and this indicates that the tidal interaction ap- 
pears weaker in the former case than the latter (because the 
curve referring to the ideal-fluid is closer to the point-mass 
curve than the polytropic curve, see below). This effect was 
already discussed in Ref. [2] and is likely due to the small 
shocks that are formed by the interaction between the outer 
layer of the stars and the external atmosphere. The polytropic 
EOS should yield a priori a more accurate evolution during 
the inspiral, when the stars are far apart, but should become 
progressively inaccurate and inconsistent when the two stars 
become closer and closer, with mass shedding and the forma- 
tion of actual shocks that are not simply due to the weak in- 
teraction with the atmosphere. For this reason we will not use 
the isentropic Q^^'s as a lower bound in our analysis, but we 
will focus only on non-isentropic evolutions, though keeping 
in mind that there is a further source of error on them. 

A natural question that comes at this stage is: what is the 
uncertainty on the determination of the Q^jii^) function that 
is due to the phase-cleaning (i.e., phase-fitting) procedure? A 
first way of addressing this issue is to measure the impact that 
changing our fiducial fitting function Eq. (28) has on Q^iiu). 
Focussing, for both models, only on our basic 1Fhr500 data, 
we computed the cleaned frequency using, besides our fidu- 
cial fitting polynomial of order n = 4 (see Eq. (28)), either a 
lower-order polynomial truncated at n = 3 or a higher-order 
one, extended up to n = 6."^ The results of these computations 
are displayed in Fig. 8 for model M2 . 9C . 12 (top panel) and 
M3 . 2C . 14 (bottom). The results are qualitatively analogous 
in the two cases. First, we see that the low polynomial order 
n ~ 3 is clearly too small, and fails to capture (when compar- 
ing it to the PN or FOB curves, which are accurate on the 
low-frequency side) the low-frequency behavior of 
By contrast, the fact that the ?i = 6 curve is barely distin- 
guishable (on the scale of the figure) from the ?i = 4 one is an 
indication of a sort of "convergence" of our fitting procedure 
as the number of a;" powers is increased. We can therefore use 
the difference between Q"^^{uj) and (5^^'*(w) as an estimate 
of the uncertainty SQ^j (w) entailed by the cleaning procedure. 



' Note that n = 5 is not meaningful as the corresponding ps term is exactly 
degenerate with (pQ. Furthermore, the use of Inx does not help, as the 
con'esponding term is nearly degenerate with 0o . 
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FIG. 9. Sensitivity of Quj to the choice of the fitting time- 
interval It for M3.2C.14. Our preferred cleaning time-interval 
/t/A/ = [1000, 2290] (central dashed-line) is compared to It /A/ = 
[1000, 2250] (solid-line) and /t/A/ = [1000, 2330] (dash-dotted 
line). See text for details. 



Computing this difference, we find that it remains of order 
unity all over the fitting frequency interval 1^. In conclusion, 
we estimate the uncertainty associated with the choice of the 
order of the fitting polynomial to be « ±0.5. Note that 
this error is rather small compared to the various numerical er- 
rors on Qi^iyj) displayed in Fig. 7, but it is only a lower bound 
on the uncertainty level S^^'^'^^Q^ linked to the cleaning pro- 
cedure. 

In particular, another relevant source of uncertainty on 
is the choice of the fitting time interval It- In Appendix B 
we explicitly discuss some rules of thumb that we follow to 
select It such that the cleaning procedure is reliable and ro- 
bust. To complete the discussion of Appendix B, we investi- 
gate (for model M3 . 2 C . 14) the modifications in brought 
by changes in the choice of It. More precisely, we modified 
the right-end point t2 of our preferred cleaning time interval 
ItjM = [/i, /2] = [1000, 2290] (see Table III) by ±40 (with 
fixed polynomial order n = 4 ). The three curves corre- 
sponding to ts = {2250, 2290, 2330} are displayed in Fig. 9. 
When comparing the cases = {2250, 2290}, we find that 
the absolute value of the difference in stays < 1 all over 
the time-interval ItjM = [1000, 1951] (corresponding to a 
frequency interval MI^^ = [0.041, 0.056]), but then grows up 
to values of order 30 near t2 = 2250. On the other hand, 
when comparing the cases = {2290, 2330} we find that 
the absolute value of the difference in stays of order 3 all 
over It- This further analysis suggests that the cleaning pro- 
cedure allows us to determine within an uncertainty level 
^cioanQ^ « 1 during most of the inspiral, with a possible in- 
creased uncertainty level w 3 near the end of the inspiral. Note 
that these levels are significantly smaller than the changes in 
the analytical Q^^'s associated to a variation of the NNLO pa- 
rameter (52 between and 100 (see next Section). 



VI. COMPARISON OF ANALYTICAL AND 
NUIMERICAL-RELATIVITY RESULTS 

A. Characterizing tidal effects from NR simulations 

Before proceeding with the NR/AR comparison it is useful 
to discuss a procedure by means of which it is possible to ef- 
fectively subtract the tidal interaction from the NR curves 
obtained so far. This procedure will then allow us to obtain 
a phase diagnostic Qj^ that, within some approximation, rep- 
resents a non-tidally interacting binary, namely a binary of 
two point-particles. As pointed out in Ref. [21], the binding 
energy of a binary system Ei,{Vl) is approximately linear in 
and it is therefore possible to subtract the tidal effects by 
combining different sets of binding-energy curves coming out 
of NR calculations. In particular, Ref. [21] computed several 
"tidal-free" binding energy curves (one curve for each com- 
bination of two different data sets) that were compared with 
the corresponding point-mass curve computed within the FOB 
approach or within non-resummed PN theory. This procedure 
allowed for both the identification (and thus subtraction) of 
systematic uncertainties in the NR data and the discovery of 
higher-order tidal amplification effects. 

Here we will generalize the approach introduced in 
Ref. [21] to the curve. In particular we assume that the 
function Q^{uj) is approximately linear in the (leading) tidal 
parameter kJ, at least during part of the inspiral, say up to 
some maximum frequency cjmax (we will use cjniax ~ 0.07). 
As a result of this assumption, we can approximately write 
Qui{uj), for each binary, as 

Q^{uj; I) = Ql{uj) + (4), Ql{u) + O ((4)^) , (35) 

where / is an index labelling some binary system. As a result, 
given the diagnostics of two different binaries with labels 
(/, J), we can estimate the two separate functions and 

Qlioj) as 



(w; J) - J) 



(36) 
(37) 



From the decomposition (35), we see that, by definition, the 
function Qj^ denotes the diagnostic of two non-tidally in- 
teracting NSs, namely of two point-like (relativistic) masses 
(and also two BHs [18, 19]). Hence, the function is 
seen to represent, within the present approximation, the effect 
of the tidal interaction on the Q^^ function. The calculation of 
both functions contains therefore important information about 
the analytical representation of tidally-interacting binary sys- 
tems. In the following we will only discuss the computation 
of the tidal-free part [uj), leaving a discussion of the prop- 
erties of QJ(oj) to a future work. 

This subtraction procedure for computing (w) can be 
first tested by using the FOB metric waveforms computed 
from binaries with compactnesses C = 0.12 and C = 0.14. 
The result of the subtraction is displayed in Fig. 10, where 
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FIG. 10. Subtraction of tidal effects: shown as a solid line is the 
point-mass BOB curve, while shown as a dashed line is the QZ curve 
obtained by inserting in Eq. (36) the tidally-modified EOB Qui curves 
shown in Fig. 3. 



we compare the point-mass (i.e., BBH) EOB curve (solid 
line) to the Q° curve (dashed line), obtained by inserting in 
Eq. (36) the C = 0.12 and C = 0.14 data of Fig. 3. The fact 
that the curves are barely distinguishable up to Mlu = 0.07 
(where the difference is AQ^^ w 1) gives us confidence that 
the procedure will be effective also with actual NR data. This 
will indeed be shown in the next Section. 



B. Inspiral: subtracting tidal effects from NR data 

We start our NR/AR comparison by computing the Q° 
function as defined by Eq. (36) from NR data using our two 
models M2 . 9C . 12 and M3 . 2C . 14 as the binaries labelled 
Ig and J in that equation. For all the comparisons carried out 
here we have limited ourselves to using the curvature wave- 
forms, although similar results can be obtained from the cor- 
responding metric waveforms. 

The results are shown in Fig. 11, which reports four dif- 
ferent Quj curves: the two tidally-modified NR Q^^ curves for 
the binaries M2 . 9C . 12 and M3 . 2C . 14 (with the asterisks 
and triangles highlighting a sample of the data on the com- 
mon frequency window), the subtracted NR (5° curve (with 
empty circles), and the point-mass-EOB (as a solid line). 
This figure illustrates at once several of the main results of this 
paper. First of all, it highlights the excellent agreement be- 
tween the cleaned NR and the analytical EOB one (cf. the 
red solid curve and the empty circles). This gives evidence 
both for the validity of the EOB description and for the ro- 
bustness of our cleaning procedure. When we compute the 
relative phase difference over the common frequency inter- 
val [0.042,0.055], we obtain the remarkably small value of 
^^EOBNR ^ ^EOB _ ^NR ^ _q 93 rad, which translates 
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FIG. 11. Subtraction of tidal effects from numerical relativity 
(curvature) Qui curves according to Eq. (36). Note the excellent 
agreement with the point-mass EOB curve in the frequency win- 
dow where M2 . 9C . 12 and M3 . 2C . 14 data overlap. The relative 
EOB-NR phase difference accumulated over this overlap interval is 
A0^°^N^ = -0.03 rad. 



into a relative difference A^^f ^^nr/^eobnr ^ o.02% 1°. 
Second, it confirms, independently of our EOB-based check 
(cf. Fig. 10), that the NR tidal effects are approximately linear 
in K2 least in the early part" of the waveform, and thus that 
they can be efficiently subtracted. Third, it illustrates the fact 
that the tidal interaction between the two objects is important 
already in the early-inspiral part of the waveform, since both 
the M2 . 9C . 12 and M3 . 2C . 14 curves are significandy dis- 
placed (by AQui ~ 10) with respect to the point-mass one. 
Fourth, such a good agreement with the point-mass EOB ana- 
lytical model (which was tuned so as to accurately reproduce 
the equal-mass BBHs) yields an independent check of the con- 
sistency and accuracy of our numerical simulations. Finally, 
we note that in Ref. [21] the procedure of subtraction, applied 
there to the NR binding energy, was giving a curve slightly 
displaced with respect to the point-mass EOB (or PN) curve. 
This displacement was interpreted as evidence of systematic 
errors in the NR simulation and prompted the introduction of 
a "correcting" procedure, which however is not necessary for 



To cross-check the consistency of both the recovery procedure of /122 from 
ip^^ and the cleaning of the phase, we carried out the same calculation also 
for the metric waveforms, finding a difference A4>^'^^^^ = +0.05 rad, 
which is consistent with the estimated error-bar A(f) = ±0.02 rad on the 
EOBNR point-mass wavefonn during inspiral [27]. 

In the following, we will refer to the frequency domain Mui < 0.06 as 
the "early-inspiral". Note that for a fiducial 1.4 — 1.4 Mq system 
Mui = 0.06 corresponds to few = 690 Hz. Note also that in the case, 
for instance, of our C = 0.14 system the frequency Mlu = 0.06 is reached 
at time t ~ 2000M, i.e., only about 5 GW cycles before merger. 
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FIG. 12. Comparison of the BOB Qi^ curves for different choices of 
the effective tidal amplification factor y4^"*'''(if) = 1 + oliu + 6l2V? , 
with the corresponding NR ones (dashed lines with open circles) 
for the two binaries considered. The dotted line corresponds to the 
"tidal free" (or "point-mass") BOB, namely, when ignoring tidal 
effects. The figure also includes the tidal-free Taylor-T4 model. 
The good visual agreement between the analytic and the numerical 
curves for 0.2 = 100 provides evidence of the need for large NNLO 
tidal corrections. The corresponding phase differences A<j)^^ — 

,EOB _ , NR ^^yg jy^ 



the present NR data. 



C. Early inspiral: evidence for large NNLO tidal effects 

We continue our analysis by focussing on the influence of 
LO tidal effects on the early-frequency part of the curves. 
We already know from Fig. 1 1 that tidal effects are impor- 
tant in such early-frequency part of the simulations, since we 
found a significant difference (of order 10) between the point- 
mass curve and the NR ones. Can these differences be ac- 
counted just by the LO tidal effects? Figure 12 shows quite 
clearly that this is not the case and that the LO description is 
not sufficient to match the corresponding NR curves (dashed 
line with open circles). Note that this is the case for both the 
M2 . 9 C . 1 2 (upper panel) and the M3 . 2 C . 1 4 binaries (lower 
panel). The difference with NR data (on the frequency interval 
[0.043, 0.057] where the data of M2 . 9C . 12 and M3 . 2C . 14 
overlap) is quantified in the first line of Table IV and is rather 
large, namely several radians. 

We next turn to analyzing the effect of NLO and NNLO 
tidal interactions. Here, we will regroup under the label of 
NLO both IPN and 1.5PN effects. As seen in Fig. 12, the 
inclusion of the NLO tidal effects (ai 1.25 [21], IPN tidal- 
radiation effects [29], and 1.5PN tail effects) has only a barely 
noticeable effect on the Q^j curve. This clearly indicates the 
need for large NNLO (2PN and higher) tidal effects, which 
we chose to parameterize by means of the effective parame- 
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introduced in Eq. (14). We then found that 
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FIG. 13. Magnitude of NNLO tidal effects: span of BOB Quj curves 
(red) with varying Q2 so as to be compatible with the various (nu- 
merical) Qu, curves (blaclc). 



choosing 0.2 — 100 yields a good match between the NR 
and FOB curves (solid line, FOB^^'^*-'), especially for 
the M3 . 2C . 14 model, for which the analytical curve is on 
top of the NR data. See also Table IV for the correspond- 
ing phase differences. The Table also indicates that if we use 
a2 ~ 130, as we did in Ref. [1], the accumulated dephasing 
on the common frequency interval [0.043, 0.057] is further re- 
duced to a fraction of a radian for both models. Note that 
the implementation of the FOB waveform and radiation reac- 
tion that we use here is slightly different with respect to the 
one of [1], which was based on Ref. [21] and thus did not 
incorporate the waveform IPN corrections [29] nor the tail ef- 
fects. This explains why we were quoting different phase dif- 
ferences (A/(/)^'-'^^'^ sa 0.1 rad) over the same interval when 
referring to a2 = 130 in [1]. However, we prefer here the 
smaller value a2 = 100 because the corresponding curve 
is, on average, closer to the NR one on the larger frequency in- 
terval [0.041, 0.068] on which we succeeded to clean the NR 
phase. 

It is important to recall that various numerical errors affect 
the computation of the NR curves, and thereby affect the 
quantitative determination of the effective NNLO parameter 
a2- For example, we have seen that the resolution extrap- 
olation (which seemed to be the dominant source of uncer- 
tainty) has the practical effect of pushing the numerical Q^^ 
curve upwards. This suggests that the value 0.2 — 100 ob- 
tained from using finite-resolution NR data is probably too 
large. To have a rough idea of the error range on a2, we com- 
pare in Fig. 13 various NR and FOB curves. More precisely, 
this figure shows two numerical black curves: a solid 
one, derived from our fiducial highest-resolution and largest- 
extraction-radius IFhr 500, and a dashed one, derived from 
the resolution-extrapolated NR data. Also reported in Fig. 13 
are three analytical curves (red lines): namely the FOB pre- 
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TABLE IV. Measuring the phase difference between NR (curvature) 
waveforms and analytic ones (from both EOB and Taylor T4 mod- 
els). The phase differences are computed on the frequency interval 
[0.043, 0.057] common to both Q^i numerical curves. From left to 
right, the columns report: the type of analytical model, the magnitude 
of the effective parameters yielding NNLO tidal corrections; and the 
dephasings A(^^^ = - <^^'^ (with X being either EOB or T4) 
for both M2.9C.12andM3.2C.14 data obtained by direct integra- 
tion of the corresponding Qui's of Figs. 12 and 14 over the common 
interval [0.043,0.057]. 



Model 


NNLO 








parameters 


[rad] 


[rad] 


EOB^° 


0.2 = 





5.04 


1.74 


gQgNLO 


0:2 = 





4.62 


1.58 


gQgNNLO 


Q2 = 


100 


1.06 


0.17 


gQgNNLO 


02 = 


130 


0.056 


-0.25 


X4LO 


T4 
012 = 


= 


6.64 


2.33 


-p^NLO 


0.2 - 


= 


6.42 


2.25 


j^NNLO 


O2 - 


= 350 


1.53 


0.15 



dictions for the three values a2 = 0, 40, 100, respectively. 
Clearly, the resolution-extrapolated curve is close to the 
analytical curve corresponding to the value 02 ~ 40, which is 
more than twice smaller than the value a2 — 100 suggested 
by our fiducial, highest-resolution NR data. It is interesting to 
note that the value 02 — 40 agrees with the preferred value 
of a2 (when using ai = 1.25) found in [21], the work that 
pinpointed the first evidence for the need of large NNLO ef- 
fects. Let us also note that, independently of the precise value 
of a2. Fig. 13 clearly shows the need for large NNLO effects, 
namely 0.2 > 40. 

Let us also recall that the other sources of numerical error 
act in various directions. For instance, non-isentropic effects 
actually act so as to effectively reduce the magnitude of the 
tidal interaction'', while the extrapolation to infinite extrac- 
tion radius acts in the opposite direction, namely effectively 
increasing the magnitude of the tidal interaction. 

In view of our incomplete knowledge of all the sources of 
error intervening in our NR waveforms, we can only conclude 
that a2 probably lies in the range 40 < 02 ^ 130, with the 
understanding that the lower values {a2 ~ 40) are preferred 
because of the expected importance of the truncation error in 
the numerical simulations. More numerical simulations with 
a more detailed estimate of the numerical error budget will be 
needed in the future to reduce this eiTor range on 02- 

Let us conclude this Section by briefly discussing the com- 
parison between the NR Qi^ diagnostics with those obtained 



Indeed the non-isentropic curve lies above the isentropic one. This is 
certainly a source of error during the early-inspiral, where the isentropic de- 
scription is a priori more accurate but some energy is channelled by shocks 
due to the interaction with the atmosphere. 




0.(M2 0.044 0.046 O.MS 0.05 0.052 0.0.54 0.056 0.058 0.06 

AIuj [curvature] 



FIG. 14. Comparison of the Taylor-T4 Quj curves for different 
choices of the effective tidal amplification factor a*"*'''(M) = 1 + 
Qi X + a2 X , with the corresponding NR ones (dashed lines with 
open circles) for the two binaries considered. The dotted line corre- 
sponds to the "tidal free" (or "point-mass") T4, namely, when ignor- 
ing tidal effects. Note that the value a J* = 350 of the dimension- 
less NNLO effective tidal coiTection parameter that best matches the 
(M3 . 2C . 14) NR data is considerably larger than in the EOB case. 
The coiTesponding phase differences A(fi^^ — <j!>^* — (j)^^' are listed 
in Table IV. 



using several versions of the Taylor- T4 approximant. IVIore 
precisely. Fig. 14 displays the following Q^^ curves: the tidal- 
free T4 model {T^^, dotted line), the LO Taylor-T4 model 
(dashed-line), the NLO (i.e., IPN) one (dash-dotted line), and 
finally the effective NNLO one (solid line), as introduced in 
Sec. IIIC above. Let us recall that the NNLO model contains 
an effective 2PN parameter, called aj^, which is a rough T4 
analog of the NNLO EOB parameter a2 and which enters the 
T4 tidal amplification factor Eq. (23). Similarly to the EOB 
case, one finds that a suitably large value of the effective 2PN 
tidal parameter aj^ can provide curves that are close to the 
numerical ones. The integrated dephasings (fp-'^ — 0^'-'^ cor- 
responding to Fig. 14 are listed in Table IV. 

A few comments are worth making on the comparison be- 
tween the EOB and T4 results. Let us first recall that, in the 
BBH case, it has been shown that the EOB description is defi- 
nitely more accurate than the Taylor- T4 one, especially when 
considering unequal mass ratios [46] or spin effects [61]. 
However, as we are considering here an equal-mass case and 
frequencies that are smaller (when considering the dimension- 
less frequencies Moj) than in the BBH case, the tidal-free T4 
phasing is quite close to the EOB one (see Fig. 12). Con- 
cerning tidal-extended models, we see that both EOB and T4 
approximations highlight the need for large higher-order tidal- 
amplification factors. When choosing one such amplification 
factor for both BNS systems (say a2 = 100 for EOB and 
aj'' — 350 for T4), a close look at the comparison of the cor- 
responding curves suggests that the EOB-predicted curves 
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FIG. 15. Comparison between BOB and NR phasing for the M2 . 9C . 12 (left panels) and M3 . 2C . 14 (right panels) binaries. The top 
panels show the real parts of the BOB and NR h22 waveforms, the middle panels display the corresponding phase differences A(/)^''''^'^^ — 
0EOB _ (f>^^^ both metric (solid line) and curvature (dashed line) waveforms, the bottom panels compare the BOB (dashed line) and NR 
(solid line) instantaneous GW frequency. The NNLO coiTections to the radial potential are carried out with the parameter a2 = 100. Note the 
agreement reached with the numerical waveform almost up to the time of the merger as defined in terms of the maximum of the GW amplitude 
(vertical dashed line) or of the contact position (dot-dashed line; see the text for details). 



are somewhat closer than the T4-predicted one to the NR 
curves. However, this, by itself, would only be a weak indica- 
tion that EOB gives a better representation of our fiducial NR 
data, especially in view of the large uncertainties discussed 
above on the actual value of the Qui'-^) functions. On the 
other hand, we consider that the need of a much larger tidal- 
amplification factor in the T4 case is an indication that the 
analytical modelling of (LO, NLO and NNLO) tidal effects 
within the EOB-resummed framework might be more robust 
than the corresponding one based on Taylor-expanded approx- 
imants. Indeed, in both cases the parametrization of NNLO 
effects involves multiplying tidal effects by a factor having a 

... ^ ^ itidal(EOB). s , , _(^) , -U) 2 

Similar structure: [u) = 1 + a\ u + a2 u ver- 

sus a^f^^{u) = 1 + aj^x + aj^x^. In addition, the quanti- 
ties u and X are numerically close to each other (both being 
close to jc?). At the end of the inspiral, 

A/w reaches numerical values of order 0.1 (i.e., 1154 Hz for 
a fiducial BNS system), corresponding to w ~ x ~ 0.136. 
For such a value one sees that the EOB amplification fac- 
tor (with OL-i = 100) remains relatively moderate'^, namely 



For 02 = 40, this amplification factor becomes 
1.25m + 40u2 ~ 1 + 0.17 + 0.74 ~ 1.91 



itidal(EOB) 



^tidai(EOB)^^^ = l + 1.25u + 100u2~ 1 + 0.17 + 1.85- 3, 
while the T4 one (with a^^ = 350) is suspiciously large, 
and is completely dominated by the 2PN contribution, namely 
a^f (u) = 1 + 1.19x + 350x2 ^ _^ o.l6 + 6.47 = 7.63. 
Another way to phrase this is to notice that the large T4 value 
aj^ = 350 is such that the 2PN contribution a^'^x^ starts 
dominating the LO term at 2- = 1/V350 ~ 1/18.7, i.e., at 
large separations r ~ 18.7il/ corresponding to rather low 
frequencies Muj = 2x^1"^ = 0.025, i.e., 285 Hz for a fidu- 
cial BNS system. Furthermore, such a large value for aj* 
works well for binary M3 . 2C . 14, but less well for binary 
M2 . 9C. 12. 



Clearly, in view of the large current uncertainties on the 
Q,^ NR curve, more work is needed to confirm this provi- 
sional conclusion. In particular, more accurate NR simula- 
tions, encompassing more compactnesses and different mass 
ratios will be needed to assess the relative merits of the EOB 
versus the TayIor-T4 description of tidally interacting BNS 
systems. 
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D. EOB/NR phasing 

So far our NR/AR comparison based on the function (cj) 
has been limited to a frequency interval which did not cover 
the last octave of frequency evolution, even if, when viewed 
in the time domain, this interval covered most of the cycles of 
the inspiral. In this section we finally focus on a phasing com- 
parison in the time domain which covers the full inspiral and 
plunge phase, up to the merger of the two NSs. Our strategy 
here will not be to explore from scratch a good range of values 
of the tidal NNLO parameter 02 values, but instead to use the 
value 0.2 = 100 suggested by our previous Qt^(cj)-analysis, 
and to explore to what extent it succeeds in providing a wave- 
form which agrees with our fiducial (highest-resolution) NR 
waveform over the full inspiral. Anticipating our conclusion, 
we will find that the EOB waveform with a2 = 100 does 
closely agree (both in phase and modulus) with the NR wave- 
form essentially up to the merger 

This is shown in Fig. 15, which compares the (real part 
of the) EOB and NR metric r/122 waveforms for the case 
including NNLO effects with 0:2 = 100. The left panels 
refer to the M2 . 9C . 12 binary, while the right panels re- 
fer to the M3 . 2C . 14 one. The top panels show the real 
parts of both the EOB and NR /122 waveforms (divided by 
the symmetric mass ratio v)\ the middle panels display in- 
stead the corresponding phase differences /S.cf)^'^^^^ {t) ~ 
^EOB^^^ — for both metric (solid line) and curvature 

(dashed line) waveforms, for completeness; the bottom panels 
compare the EOB (dashed line) and NR (solid line) instanta- 
neous GW frequency. The least-squares phase alignment has 
been performed on the time interval [tijifl] /A/ = [250,3300] 
for the M2 . 9C . 12 binary and [i^, <fl]/Af = [250, 2250] for 
the M3 . 2C . 14 one. 

The two vertical lines (dot-dashed and dashed) indicate the 
"end of the inspiral phase", as defined either within the EOB 
analytical framework (dot-dashed line) or by using NR infor- 
mation (dashed line). Note that we call here simply "inspi- 
ral" what was called "insplunge" in previous EOB studies, 
namely the union of the inspiral and (when it is reached before 
merger) of the plunge. More precisely, the dashed line indi- 
cates the NR-defined "merger", i.e., the time (computed from 
the NR data) at which the modulus of the metric waveform 
reaches its first maximum. On the other hand, the vertical 
dash-dotted line indicates the EOB-defined "contact" between 
the two NSs'"^. Such a formal contact moment was introduced 
in Eqs. (72) and (77) of Ref. [21], by a condition expressing 
that the EOB radial separation R becomes equal to the sum of 
the tidally deformed radii of the two NSs, namely 

j^contact ^{l^h^ e^(i?™"*^^')) Ra + {a ^ b} , (38) 

where eA{R) = MbR\/{R^Ma) is the dimensionless pa- 
rameter controlhng the (LO) strength of the tidal deformation 



of the NS labeled A by its companion B and where ' 
is the shape Love number [18, 62]. A recent study of the 
tidally induced shape deformation of BHs [62] has shown that 
the BH shape Love number was a function of the sepa- 
ration r, which increased as r decreased (and u increased). 
This behavior is similar to the behavior of the (effective) 
quadrupole Love number k^{u) = fc2(l + otf^u + 0^2^ v?'), 

(2) (2) 

where both a\ and were found to be positive [21]. One 
would need a special study devoted to the comparison of the 
EOB-predicted NS shape deformation to NR data to investi- 
gate in detail the u dependence of the analogous h'^{u) — 

^2(1 + I'l^u + 72^-' u^). Leaving to future work such a study, 
we will here replace the u-dependent effective shape Love 
number ^^{u) by a constant, chosen such that the EOB- 
predicted contact happens before the NR-defined merger for 
the two BNS systems we consider We found that h!^ ~ 3 
works, and this is the value we will use to replace /i^ and 
/i^ in the contact condition written above'^. An important 
point to note is that our (EOB-based) analytical definition of 
contact allows one to analytically predict a complete inspiral 
waveform, including its termination just before merger 

Figure 15 shows that the agreement in the time domain be- 
tween the analytic EOB description and the fully numerical 
one is extremely good essentially up to the merger More 
precisely, the match between the two descriptions is excel- 
lent both in modulus and in phase, with a dephasing of order 
A(/) = ±0.1 rad during most of the long inspiral phase. It 
is only during the last 100 A/ before contact that the dephas- 
ing grows significantly. One should note that this excellent 
EOB/NR agreement holds for both binaries M3 . 2C . 14 and 
M2 . 9C . 12, and has been obtained by tuning a single tidal- 
amplification parameter 

Clearly the results presented here give only a first cut at 
these issues. More NR/AR comparisons are needed to con- 
firm our findings and to determine the most effective value 
of 0.2- With sufficiently accurate NR data one can hope 
to determine not only the effective tidal-amplification factor 
^tidai^y-j _ 2^ ^ _|_ but also the precise separation- 
dependence of A^l'^^'^{u). This would allow one to extend the 
EOB description right up to the merger. 

VII. CONCLUSIONS 

We have presented the first comprehensive NR/AR com- 
parison of the gravitational waveforms emitted during the in- 
spiral of relativistic BNSs as computed via state-of-the-art 
numerical-relativity simulations and as modelled via state-of- 
the-art analytical approaches. Overall, the work reported here 
and our findings can be summarized as follows. 

1 . We have considered the longest to date numerical simu- 
lations of inspiralUng and coalescing equal-mass BNSs 



Note that the styles of the corresponding merger and contact vertical lines 
as depicted in the two panels of Fig. 2 of Ref. [1] are inverted with respect 
to the text there. See the arXiv version for the conect figures. 



A similar approach was taken in [1, 18], with a less conservative value 
h"^ = 1. Let us recall that the computation of the infinite-separation shape 
Love number h2 = h2^[u = 0) of NSs gives values of order unity [18]. 
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modeled either with an ideal-fluid or a poly tropic EOS. 
Because tidal effects are most sensitive to the stellar 
compactness, we have considered two binaries with ei- 
ther a small compactness of C = 0.1199 or with a large 
compactness of C = 0.1396. The parts of the wave- 
forms relative to the inspiral cover between 20 and 22 
cycles and have been studied to isolate possible sources 
of error, such as non-isentropic evolutions, finite-radii 
GW extractions, and the use of finite resolutions. For 
the model with the highest compactness, the first two 
sources of errors lead to a total error-bar in the GW 
phase of ~ ±0.15 rad, while the error coming from 
a finite resolution indicates an accumulated phase error 
of A0c:i ±0.54rad. 

2. We have used the function Qi^{uj) = lo'^ /6j as a use- 
ful diagnostic of the physics driving the evolution of 
the GW frequency uj. The calculation of this quan- 
tity is however challenging when made from the early- 
inspiral part of the NR waveforms, as the latter is af- 
fected by a series of contaminating errors. We have 
filtered out these errors by fitting the NR phase evolu- 
tion <l){t) with a simple analytical expression that re- 
produces at lower order the behavior expected from 
the PN approximation. We have compared the various 
Qtj's obtained from different data to estimate the er- 
ror range entailed by comparing analytical predictions 
to our highest-resolution, largest-extraction-radius NR 
data. 

3. Using the estimated Q^i'-^) function we have shown 
that it is possible, at least for frequencies Mijj < 0.06 
(i.e., /gw < 700 Hz for a fiducial 1.4 A/© - 1.4 M© 
BNS system), to subtract the tidal-effect contribution 
from the NR waveforms and consistently match this 
with the expected EOB model for point particles which 
has been successfully matched to BBH simulations. 
The ability to perform this match accurately provides 
us with an independent validation of the quality of our 
numerical results as well as with a confirmation that the 
function Qui{uj) is approximately linear in the (leading) 
tidal parameter k^. 

4. The comparison of analytical predictions with NR data 
shows that tidal effects are significantly amplified by 
higher-order (NNLO) relativistic corrections even in the 
early inspiral phase. These NNLO tidal corrections are 
parameterized within the EOB approach by a unique 
(effective, 2PN) tidal parameter a2- Although the most 
precise available at the moment, the quality of the NR 
data is such that we can only constrain the actual value 
of (52 to be in the range 40 < a2 ^ 130. 

5. Once a single choice for a-z is made, the EOB-predicted 
waveforms agree (both in phase and in modulus) with 
the NR ones (for both BNS systems) within their error 
bar and essentially up to the merger 

6. Finally, we have also compared the NR phasing with the 
one predicted by a non-resummed Taylor- T4 PN expan- 



sion, completed by additional tidal terms. If one uses 
only the currently known analytic T4 tidal terms, the T4 
model dephases (when C = 0.12) by more than 27r rad 
already at the GW frequency Muj = 0.057, which is 
about twice smaller than the GW frequency at merger 
(we recall that Muj = 0.057 corresponds to 658 Hz 
for a fiducial 1.4 Mq - 1.4 M© system). On the other 
hand, a good match (for both compactnesses) with the 
NR phasing is possible if one allows for a T4 analog of 
the EOB Q!2 parameter, i.e., an (effective) 2PN ampli- 
fication of tidal effects. The corresponding parameter 
aj"' ~ 350 is suspiciously large, works well for binary 
M3 . 2C . 14 but less wefl for binary M2 . 9C . 12, and 
dominates the amplification of tidal effects already at 
frequencies Mlu ~ 0.025 (corresponding to 285 Hz). 
This seems to suggest that the EOB-based representa- 
tion of tidal effects is more reliable than the Taylor- T4 
one. 

In summary, the work presented here opens new avenues 
to the important synergy between numerical and analytic de- 
scriptions of inspiralling compact-object binaries in general 
relativity. For the first time we have shown that an analytic 
modelling is possible also for objects which cannot be treated 
as point-particles and for which, therefore, tidal effects rep- 
resent important corrections. Although the results presented 
here are very encouraging, a number of improvements are 
needed on both the numerical and the analytical sides. On 
the numerical side, higher resolutions and better measures of 
the convergence rates (which are particularly challenging in 
non-vacuum simulations) are needed to decrease the numeri- 
cal phase errors to and reach firm conclusions about the tidal 
contributions to the phasing. On the analytical side, higher- 
order PN calculations are needed to better determine the form 
of the NNLO corrections. Both of these goals will be the sub- 
ject of our future work. Hopefully, progress on both fronts 
will enable us to determine the crucial tidal-induced dephas- 
ing function A'"^^V(<^) with an accuracy sufficiently high to 
extract reliable information on the EOS of matter at nuclear 
densities 
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Appendix A: Computing metric waveforms from 4>4 

We discuss here the details of how to accurately derive 
the metric waveforms /i+,x from the numerically computed 
curvature waveforms We first recall that the procedure 
outlined in Ref. [46] consisted essentially of three steps, (i) 
First one performs the double integration of ^l"' starting at 
t = with the integration constants set to zero; this amounts 
to defining 

h'j^it)= f dt'ii^{t'), (Al) 
Jo 

/i^™(t)= ( dt'hi"'{t'). (A2) 
Jo 

The provisional metric waveform /ig™ differs from the "exact" 
metric waveform (6) (integrated from past infinity) by a linear 
function of t, say 

/i^Q™ (t) = h'-^ (t) + «exactt + /3oxact • (A3) 

(ii) The second step consists in obtaining an estimate of the 
two (complex) integration constants (ctoxact, /^oxact) that en- 
ter the metric waveform (6) by fitting over the full simula- 
tion time interval the [t > Q)-integrated waveform (A2) to 
a linear function of t, say /ijj""^* = at + j3, where a and 
/3 are complex quantities, (iii) The third and final step of 
the procedure of Ref. [46] consisted in subtracting the lin- 
ear function at + (3 from Hq™ so as to define an approx- 
imation to the {t > — o3)-integrated metric waveform, say 

Here we will use a "new" (three-step) procedure, which 
starts with the same step (i), but modifies both steps (ii) and 

(iii) so as to get a better approximation to the exact metric 
waveform. First of all, we define an "adiabatic-like" approxi- 
mation to the metric waveform, 

h,^{t) ^ -^f^ (A4) 

and use this to define 

h'j-{t) = h'j-{t)^h„n{t). (A5) 

As h^''"{t) is approximately equal to h^"^{t) (because of the 
approximately adiabatic nature of the inspiral), we see from 

Eq. (A3)that/l^'"(t) = /lfm(t)-/lfm(i)+Q!exacti+/3exact will 

be closer to the unknown linear function Oexacti + /3exact than 
/iQ™(t) was. Therefore, the next step is to perform the linear 
fit on this /lo™ rather than on hf^'^{t) itself. Then, the last 
step (iii) consists, as in the past, in subtracting the resulting 
improved linear fit at + f3 from the {t > 0)-integrated metric 
waveform h^^'"-{t). 



In addition, let us note that we perform the fit not on the 
whole time interval, but rather on a restricted time interval 
that cuts away the first cycles of the waveform. Finally, after 
doing several tests, we realized that the entire procedure leads 
to a physically more reliable metric waveform (see below) if 
/iQ™(t) is fitted not to a simple linear function, but rather to a 
quadratic'^ one, hl'''"^-^\t) = -ft'^ + at + l3. 

As emphasized in Ref. [46], we accept the integrated wave- 
form if and only if its modulus exhibits a monotonic growth 
in time during the inspiral, consistently with the expected cir- 
cularly polarized behavior of the metric waveform (as well 
as the curvature one)"^. Figure 2 displays the metric wave- 
forms (both for the M2 . 9C . 12 (left) and the M3 . 2C . 14 
(right) models) obtained using this improved procedure. The 
time intervals where we fit the waveforms to get 

^quad-fit 

Start respectively at ti/M 294 (model M2 . 9C . 12) and at 
ty/M = 677 (model M3 . 2C . 14). Note how the modulus of 
both models exhibits a smooth monotonic behavior in time. 



Appendix B: Cleaning tlie GW phase and Qui curves 

We next provide more detailed information about the clean- 
ing procedure of the NR GW phase advocated in Sec. IV and 
used to drive NR/AR comparisons. As we said in the main 
text, the final goal is to fit away the high-frequencies oscil- 
lations in the GW phase (j> so as to get a clean and smooth 
Quj curve, Eq. (24). We recall that the idea is to fit with 
an analytic expression that is modeled on the PN expansion. 
Defining the quantity 

riy -1-1/8 
x(i,^c) = {g(tc-t)} , (Bl) 

one then fits the NR phase with an expression of the form 
9 

(t)=--x-'' (l+P2X^ +P3X^ +P4x'^ + ...) +00, (B2) 

where t^, <\>o, and the p/s are free coefficients to be deter- 
mined by the fit. Note that tc can be thought of as defining 
a formal "coalescence" time. There are two deUcate (corre- 
lated) points: (i) how many powers of x [possibly including 
also ln(a;) terms] one has to include in Eq. (B2), and (ii) 
on which (time) interval ItjM = (ii, ^2) the approximate de- 
scription of given by Eq. (B2) (and consequently of Q^) is 
reliable. The procedure to select the "best" time interval and 



We think that such a quadratic fit is needed for absorbing several effects that 
"pollute" the waveform, notably finite-extraction-radius effects, remnant 
junk radiation, etc. In this respect, we also mention that Ref. [45], in the 
context of non-spherical star oscillations, found that a quadratic polynomial 
used in the recovery of /120 from il)^ was a necessary choice to find a good 
agreement both with Abrahams-Price metric extraction and perturbative 
waveforms. 

Note however that small-amplitude, high frequency "ripples" are still 
present in the modulus. Their origin is however essentially numerical, as 
they are also present in the modulus of . 
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FIG. 16. Testing the fit of the GW phase of the M2 . 9C . 1 2 simulation. The top-left panel shows the time evolution of the frequency, computed 
from the metric and curvature waveforms. The bottom-left panel shows the deviation of the cleaned phase evolution with respect to the raw 
data; note that they average to zero. The right panels show the comparison of the frequency evolution of the cleaned and raw wavefoms, for 
the curvature (top) and metric (bottom) waveforms. 



to consistently assess the quality of our cleaned curves can be 
summarized as follows; 

1 . The initial time ti is chosen so as to eliminate as much 
as possible the most noisy part of the curvature fre- 
quency. In practical terms, this meant cutting at ti = 
1200 for M2 . 9C . 12 data and h = 1000 for M3.2C.14 
data. This is illustrated in the top-left panels of Fig. 16 
(for M2 . 9C . 12 data) and of Fig. 17 (for M3 . 2C . 14 
data), which show the curvature (dashed line) and met- 
ric (solid line) instantaneous GW frequency oj. In both 
plots, the first vertical line identifies the location of ti. 

2. For a given order of the polynomial, we found the right 
end, t2, of the time window essentially, by trial and er- 
ror, monitoring the behavior of several quantities. In 
particular, (i) we checked that the cleaned oj visually 
"averages" the raw lu, for both -04^ and /122 data. This 



is illustrated in the top-right and bottom-right panels of 
Figs. 16-17, the raw data appearing as dashed lines, the 
cleaned data as solid lines. Then, (ii), we require that 
the phase difference ^^ican _ ^Raw averages to zero, 
which indicates that we have subtracted all the "secu- 
lar" trends by means of our polynomial fit. The quan- 
tity A^cicanRaw ^ ^cioan _ ^Raw (^^^^ curvature and 

metric) is displayed in the bottom-left panel of Figs. 16- 
17. The fact that it averages to zero is the indication that 
our fit caught the "secular" behavior of the phase, aver- 
aging away both (numerical) low-frequency and high- 
frequency oscillations. 

3. For a fixed time window, the inspection of /^^CieanRaw 
is also crucial for choosing the order of the polynomial 
in X, which we set to be of fourth-order A 3rd-order 
one is clearly not enough to get the right trend of the 
frequency (and thus of Q^) up to the end of our pre- 
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f erred interval. 

4. To better select the end t2 of the time window, we found 
it useful to monitor the difference between the curva- 
ture and metric Q^'s, namely AQ^"™ = Q^^'vaturc _ 
Qmetnc ^y^e typically choose the value of tj^ in such a 
way that AQJ^^'" is always smaller than 0.2 on the fre- 
quency interval corresponding to It/M. This value can 
be estimated by comparing curvature and metric Q^j's 
within the EOB. For example, for the NNLO model 
with (52 = 100 one checks that AQ^"™ < 0.2 when 



uj e [0.035,0.055] forC = 0.12, and AQ^-"' < 0.2 
when UJ e [0.035, 0.063] for C = 0.14. This gives us 
an idea of the level of AQ"^"™ that we can accept from 
our cleaned NR curves, so that we can choose the fitting 
time window accordingly. 

In conclusion, to obtain the central NR-cleaned Q^^ curves 
labelled IFhr500 used in the core of the paper, we fixed 
tji/M = 3366 for the M2 . 9C . 12 phase and tn/M = 2290 
for the M3 . 2C . 14 one. The time intervals (and the corre- 
sponding frequency ones) used to clean the other NR phases 
are also listed in Table III. 
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